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PART I: AN OVERYIEV 


1. Introduction 

At present, metal matrix composites are the primary materials for use in 
the propulsion systems of space vehicles due to the light weight/high perfor- 
mance requirements. During its service life, the composite will experience, 
to a large extent, severe thermal and mechanical loading cycles under high 
temperature environments. As a result, such composite components undergo 
considerable inelastic deformations, leading to phenomena like ratcheting, and 
eventually to their fatigue failure. Obviously, effective utilizations of 
high temperature composites require proper design technology, including life 
prediction methodologies, for "typical" composite structural components; such 
as turbine blades and combustor liners in the form of composite laminates with 
flat (plate- type) or curved (shell- type) geometries. 

A valid assessment of the structural integrity, reliability, and life 
expectancy of these components requires the development of an improved numeri- 
cal capability for their complete, "global- local" (also called progressive- 
failure) analysis. This in turn requires consideration of (a) generalized 
material behavior, and (b) analytical solution procedure. It is noted here 
that for such a development, a general framework of wide applicability is 
still lacking. The research work reported here is concerned with several con- 
tributions to this end. 

The term "generalized material behavior" refers to the multiaxial consti- 
tutive equations which describe the basic characteristics, and experimentally- 
observed phenomena, of the composite material subjected to complex thermo- 
mechanical load cycles. These are the most fundamental relations required in 
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any analysis; in fact, their "accuracy" determines the "quality" of the entire 
numerical results. Ideally, these constitutive models should then provide 
accurate mathematical representations of the pre- failure, failure, and post- 
failure (or damaged) materials response modes. 

However, in view of the extremely complex nature of the many phenomena 
exhibited by composite materials at elevated temperatures (e.g., strong 
initial anisotropy; inelasticity; time- and temperature- dependent effects like 
creep, relaxation, and recovery mechanisms; failure in the form of fiber 
breakage or matrix- dominated ductile rupture; constituent matrix- fiber inter- 
actions such as debonding and delamination failure, etc.), a concensus favor- 
ing a particular approach for developing a "comprehensive" constitutive model 
for composite deformation- strength behavior is not yet in evidence. Neverthe- 
less, progress is currently being made on several fronts by various groups at 
NASA Lewis [e.g. 1-9]. Considering the "output" composite constitutive ideal- 
izations resulting from this research activity, an immediate and important 
issue in this connection concerns their use in nonlinear structural analysis. 

The term "analytical procedure" used above refers collectively to all 
matheraatical/numerical aspects of calculations needed to obtain a solution. 
To this end, the finite element method is often utilized as a general frame- 
work. The present work is mainly concerned with issues related to accuracy, 
reliability, and efficiency of algorithmic strategies employed in the (ma- 
terially/ geometrically) nonlinear finite element solutions for laminated 
composite plates and shells under static/dynamic loadings. 

More detailed discussions are given below on specific areas and problems 
investigated in the subsequent parts of the report. 
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2. Composite Constitutive Modeling 

Metal matrix composites are typically multi- phase materials, comprised of 
"stiff" reinforcing fibers, metal matrix (resin), as well as their inter- 
phases. The mathematical characterization of their "macroscopic" thermo- 
mechanical deformation/strength properties should therefore ideally be based 
on a consideration of the more basic "microscopic" properties of the indivi- 
dual constituents and their interactions; i.e. in the spirit of the so-called 
"mesomechanics" approach [6] . 

In one such micromechanics approach, the overall response of the aggre- 
gate composite is "derived" using the rule of mixtures, and often utilizes 
some simplifying assumptions regarding the fibers’ geometry and packing, for a 
representative volume element of the composite [1,8,9]. For example, this has 
led to the development of a set of "simplified" nonlinear constitutive rela- 
tionships by Chamis and his colleagues at NASA Lewis [e.g. 1]. 

Other more comprehensive developments of micromechanics- based inelastic 
constitutive models for composites at elevated temperatures are currently 
underway at NASA Lewis using the method of homogenization. For instance, 
based on the periodicity of the composite microstructure, a combined experi- 
mental/finite element effort is made by Onat and Leckie [7] to formulate con- 
stitutive/damage models for the "equivalent" homogenized composite material. 
In addition, homogenization techniques based on Green functions/Fourier series 
are employed by Valker [3] in developing viscoplastic models for general 
periodic/nonperiodic heterogeneous composites, which can also account for 
surface effects in thin structures. 
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However, in an alternative approach for composite material modeling, the 
microscopic effects can be "averaged" at the outset, and the phenomenological 
(experimentally- observed) aspects of the composite can then be idealized as 
for an initially- anisotropic continuum. Various classes of such continuum- 
based constitutive representations have been recently developed at NASA Lewis; 
i.e., unified viscoplastic models by Robinson and co-workers [e.g. 4, 10]. 

In view of their generality and ability to account for many important 
cyclic load/time- dependent phenomena, these multiaxial viscoplastic equations; 
e.g., as suggested in [4] for transversely- isotropic composites, are presently 
among the most well- developed and promising models for practical applications. 
They are therefore used here as a base- form for the finite element capabili- 
ty developed. Note that further extensions/refinement of this basic formula- 
tion are also currently under investigation; e.g. to include high- temperature 
damage (creep damage, low- and high- cycle fatigue damage, etc. [7]). Ve also 
note that the recent research efforts on homogenization in [3,7] will most 
likely lead to constitutive macromodels of this same basic structure. 

Vith regard to the above viscoplastic constitutive modeling, a major part 
of the present work is concerned with the study of various computational/ 
implementation aspects associated with their usage in large-scale finite 
element analysis [e.g. 5, 11-17]. More specifically, the main objectives here 
are: 

(1) Development of efficient/accurate integration schemes [e.g. 11-14] 
for the resulting system of stiff nonlinear differential equations 
using explicit/implicit methods with error control; 
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(2) Design of special procedures for the structured- coding organization 
of the material model implementation into the finite element soft- 
ware NFAP [18] , such that it becomes not only transparent but also 
immediate to constitutive researchers at NASA Lewis for future 
work. 

(3) Development of a "refined" automatic local/global time incrementing 
algorithm [e.g. 15-17] capable of handling general loading; e.g. 
monotonic or low- rate time varying load histories, as well as trans- 
ient or high- cycle thermomechanical loadings, etc. 

Detailed descriptions of the above items and results of the numerical 
tests for the schemes developed are given in Part III of this report. 

Several other aspects that are not addressed in this report, but will be 
considered in our future research, are associated with the introduction of the 
"damage" response [19] in the viscoplastic model. It is now well-known that 
this presents unique numerical difficulties and calls for a careful investiga- 
tion of the computational procedure utilized. For example, a typical con- 
tinuum-damage model is likely to eventually lead to some sort of "strain- 
softening" [20,21]. This may cause the finite element solution to exhibit 

pathological behavior in the form of oscillations and strong element- size/ 
mesh sensitivity, unless suitable "localization limiters" are introduced in 
the numerical model [21]. Also, depending on the particular continuum- damage 
theory used, the resulting "damaged" material- moduli matrix may become "un- 
symmetric" [20,22]; i.e., an efficient symmetrization procedure may thus be 
needed. 
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Finally, there are cases in which composite structures may fail due to 
excessive or "moderately- large" viscoplastic deformations (e.g. strains of the 
order of a few percent, as opposed to the infinitesimal- strain assumption 
typically used) . One crucial question in such cases is how should these 
"moderate" inelastic strains be treated in the finite element analysis? For 
example, the inelastic deformations appear as an additive term in the 
constitutive equations of the unified viscoplastic theory, thus leading to the 
so-called initial strain approach which is typically employed in finite 
element analysis. This is equivalent to an explicit or constant- stiffness 
formulation for "small- strains" viscoplasticity. It is known that the expli- 
cit method may present numerical instability or slow convergence problem, 
especially when the composites experience large strains. In this latter case, 
a need therefore arises for the development of variable- stiffness or implicit 
methods [e.g. 23], which has been traditionally used for time- and rate- 
independent plasticity. This in turn requires careful studies of the possible 
ways for deriving appropriate material "tangent- stiffness" matrices for 
viscoplastic models. 

3. Laminated Composite Plate/Shell Elements 

A very important ingredient in the nonlinear analysis of composite struc- 
tural components is the use of suitable plate/shell elements representing 
structural action of thin/moderately- thick/thick laminated composites. It is 
presently well-established that the most effective approaches for geometric 
modeling of plates and arbitrary curved shells are: (a) degenerated shear- 
flexible (middle- surface) elements for thin and relatively- thick situations; 
and (b) two- surface elements (3-D solids) for the thick regime. Despite 
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their popularity, displacement- based models of these types are also known to 
exhibit several difficulties [e.g. 24, 25, 27] in applications to isotropic as 
well as anisotropic problems (e.g., shear/membrane locking for thin struc- 
tures, kinematic deformation modes for reduced integration, etc.). 

Alternative formulation approaches have therefore been advocated, e.g. 
using various hybrid/mixed methods. In particular, a class of effective and 
simple hybrid/mixed plate and shell elements have recently emerged from our 
previous research work with NASA Lewis [25-30]. More specifically, a large 
variety of critical test problems, for both linear [25, 27] as well as geo- 
metrically- nonlinear [30] situations, have clearly demonstrated the robustness 
and accuracy of our simple quadrilateral element (HMSH5) for isotropic- 
elasticity. A major portion of the present work was thus devoted to mixed 
element- technology development applied to composites, for both static and 
vibration analyses, as detailed in Part II of the report. 

In such a development of mixed composite elements for plates and shells, 
at least three numerical problems can be anticipated. First, because of their 
anisotropy there will be strong interactions among various stress and strain 
components in a composite element, which may induce the so-called locking 
effect. This needs careful examination. For example, an attempt has been 
made by Spilker [32] to formulate a locking- free 8-noded hybrid stress element 
for laminated plates. However, in that approach, independent stress functions 
were assumed for each individual lamina. This leads to a very large number of 
independent stress (strain) parameters per element, thus rendering the analy- 
sis prohibitively expensive. Instead, we propose to use a fixed set of stress 
(or strain) functions for a "designated" number of laminae (plies or layers). 
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This "designated" number has to be determined on the basis of numerical study. 
In fact, the extreme assumption of utilizing layer- number- independent strain 
interpolation for the laminate element, is utilized here, and was shown to be 
successful in all the test cases investigated in Part II. This is certainly 
of great advantage from the standpoint of computational efficiency. 

Second, there is a need in using an efficient "through- thickness" inte- 
gration scheme to obtain the overall constitutive properties of the laminated 
composite from the corresponding (viscoplastic rate) equations for the indi- 
vidual plies (laminae). The familiar Gauss quadrature rule for isotropic 
materials is not suitable for laminated composites since it fails to capture 
the material behavior of each lamina. For numerical expediency, we developed 
an "interpretive" scheme for thickness integration that is combined with 
conventional Gauss quadrature for in- plane lamina integration. That is, a few 
sampling points in the thickness direction are preselected to calculate the 
material matrices for an equivalent unidirectional material. Then, the cor- 
responding quantities for each lamina are interpolated according to its 
thickness- location and fiber- orientation. It is believed that this scheme 
combines the attributes of computational efficiency as well as numerical 
accuracy. 

Finally, from the viewpoint of general applicability to arbitrary curved 
shells (i.e., nonflat/nonrectangular elements), the need will arise for the 
appropriate definition of "fiber- orientations" or "material- symmetry" axes for 
the "equivalent" material at each in- plane (lamina) integration point. To 
this end, we have implemented a simple "shape- function" interpolation (iso- 
parametric-type) scheme based on their values associated with the nodal points 
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defining the element geometry; i.e., the so-called nodal fiber or material 
triads. Note that for rectangular- element meshes for flat plates, this can be 
reduced to simply specifying (as input) the angle of inclination of material- 
symmetry axes with respect to the global axes of the composite plate. 

The various aspects above have been incorporated in the mixed formulation 
of the composite element HMSH5, and its applications to a number of test cases 
are described in Part II of the report. Emphasis is placed on demonstrating 
its response characteristics with regard to: 1) mesh- convergence properties; 
2) accuracy of stress and displacement prediction in static analysis; and 3) 
frequency and mode predictive capability for vibration analysis. 

4. Fully- Nonlinear Analysis Capability for Shell Applications 

In addition to material nonlinearities, geometric nonlinearities in the 
form of 1) large rotations , and 2) large or non-infinitesimal strains, must 
be also considered in the development of a fully- nonlinear solution procedure 
of plates and shells [e.g. 30, 31]. For example, item (1) is crucial in 
studies related to any buckling phenomenon (e.g. creep buckling). Signifi- 
cance of (2) appears in the analysis of failure or damage due to excessive 
deformations. In these latter cases, moderate strains (on the order of 3-57.) 
will be typically involved. Based on our previous experience [30] , these 
should be properly accounted for in the iterative solution scheme to improve 
its convergence properties. 

In this final part of the investigation, the general framework of fully- 
nonlinear curved shells is developed for the mixed elements. Adopting an 
updated- Lagrangian approach, a consistently- linearized form of the incremental 
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modified Hellinget/Reissner mixed variational principle is utilized to derive 
the element nonlinear governing equations. Emphasis is placed on devising 
effective solution procedures to deal with 

(1) element configuration- update in the presence of large rotations in 
space; 

(2) accurate integration of spatial stress and strain fields. 

A rather comprehensive account of the various aspects of this development 
has been recently compiled in a report [30], which was provided to NASA- Lewis 
earlier this year (1989). Therefore, only a summary of the essential points 
and results in some validation tests will be included in Part IV of this 
report . 

5. Objectives. First- Year Accomplishments, and Outline 

In summary, the overall objectives of the present work are given as 
follows: 

1) To develop robust and effective laminated plate/shell elements based on 
the mixed method for the analysis for static and vibration problems. 

2) To develop efficient implementation and time- stepping/numerical- integra- 
tion schemes for unified viscoplastic constitutive models for high- 
temperature composites. 

3) To develop a general framework for the fully- nonlinear structural analy- 
sis of shells. 

Considering the above three items, our first- year accomplishments are 
summarized as follows: 

Vith regard to item (1), we have completed the formulation and testing of 
an "extended" hybrid- mixed composite element HMSH5 for static and vibration 


10 



analyses. Emphasis is placed in this development on: i) robustness, ii) com- 
putational efficiency, as well as iii) general applicability to arbitrary 
curved shells, and several noteworthy aspects are included here. With regard 
to (i) a careful selection is made for the polynomial functions in the assumed 
strain field, thus leading to stiffness formulation that is kinematically 
stable, free from locking for large "thickness" aspect ratios, relatively in- 
sensitive to geometric distortions, etc. In addition, with a view toward 
(ii), and in contrast with other previously- developed hybrid composite plate 
elements, we utilized a fixed set of parameters; i.e., a layer- number- indepen- 
dent assumption, for the strain field discretization, together with an effi- 
cient "through- thickness" integration scheme for the stiffness calculations in 
terms of the anisotropic elastic/viscoplastic properties of the laminated ele- 
ment. A simplified mass matrix was also used, thus obviating the need for 
expensive "dynamic condensation" of its fifth node’s degrees- of- freedom. 
Finally, from the viewpoint of general applicability to arbitrary curved 
shells (i.e., nonflat/nonrectangular elements) in (iii), appropriate defini- 
tions of "fiber- orientations" or "material- symmetry" axes for the "equivalent" 
material at each in- plane (lamina) integration point were made through the 
implementation of a simple "shape- function" interpolation (isoparametric- type) 
scheme based on their values associated with the nodal points defining the 
element geometry; i.e., the so-called nodal fiber or material triads. 

In connection with research areas in item (2) above, the major contribu- 
tions of the first- year study consist of: i) four different viscoplastic 
constitutive models (for isotropic and anisotropic responses), and the asso- 
ciated integration methods, are coded as a separate module in NFAP so that any 


11 



future extension or coding modification can be done easily by a user, ii) "re- 
duced" forms of viscoplastic models wqre made readily available for applica- 
tions to various structural components represented by a 2/D continuum, or 
plate/shell elements, iii) further understanding of the numerical schemes 
employed for practical analysis, and iv) development of a "refined" automatic 
local time incrementing algorithm capable of handling general time- varying 
load histories. 

Finally, in the first- year investigation associated with item (3), the 
general framework of fully- nonlinear analysis of curved shells has been devel- 
oped for our mixed finite elements. Adopting an updated- Lagrangian approach, 
a consistently- linearized form of the incremental modified Hellinger/Reissner 
mixed variational principle is utilized to derive the element nonlinear 
governing equations in the form of tangent ( variable ) stiffness solution. 
Emphasis is placed on devising effective solution procedures to deal with: i) 
element configuration- update in the presence of large rotations in space; and 
ii) accurate integration of spatial stress and strain fields when large non- 
infinitesimal strains occur. Vithin this framework, the newly- developed non- 
linear element HMSH5 was tested and shown to exhibit improved convergence 
characteristics compared to other available elements. 

Following the overview given in this first part, the remainder of the 
report is conveniently divided into three separate parts, II, III and IV, des- 
cribing the research work completed in each of the above three areas (1) , (2) , 
and (3), respectively. Each part contains a review of pertinent literature, 
detailed theoretical and algorithmic developments, results of numerical test 
problems for validation, as well as conclusions summarizing the research items 
planned for future work. 
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PART II; MIXED ELEMENTS FOR LAMINATED PLATES AND SHELLS 

1. Introduction 

In order to analyze and design some of the complex composite structural 
applications, finite element techniques are now almost invariably used in 
many, if not all, phases of the design process. In addition to the design 
phases, recent interest has also begun to focus on the failure and various 
damage mechanisms associated with composite materials, i.e. matrix cracking, 
delamination/debonding, ply failure, etc. These types of applications demand 
accurate stress predictions. 

Because of the above mentioned points, considerable interest and effort 
has focused in developing accurate and efficient composite finite elements. 
In this connection, we summarize here the results of our recent research 
activity on the development of an effective quadrilateral plate/shell element, 
of the hybrid/mixed type, for the static and dynamic analysis of anisotropic 
elasticity. 

2. Background and Literature Review 
2.1 Composite Finite Elements 

A laminate is typically constructed of a series of plies stacked and 
bonded together to form the laminate. Each of the individual plies are com- 
posed of two distinct constituents, fibers and surrounding matrix material, 
which possess markedly different properties. In particular, due to the 
difference in shearing stiffness of the matrix and fibers, high ratios of 
elastic moduli to shear moduli are developed. Thus, the cross-sectional 
warping, i.e. shear deformation, of the laminate is dependent on the ply 
orientations in addition to the plate/shell thickness. Because of this, it is 
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known that shear deformation has a more pronounced effect in laminated 
structures than in the case of isotropic structures. 

In classical lamination theory (CLT) the usual Love-Kirchoff assumptions 
of plane sections remaining plane are in effect, thereby neglecting shear 
deformations totally. Obviously, this type of theory is of no value except in 
very thin plate or shell applications. A theory used in thick plates is that 
attributed to Reissner-Mindlin type of assumptions. In this case, a constant 
shear angle through- the- thickness of the plate is assumed and has been refer- 
red to as the constant shear angle theory (CST). Typically, some form of 
shear correction factor, k , is chosen to account for the through- the- thickness 
shear deformation. These types of Reissner-Mindlin elements are those typi- 
cally used for composite analysis. 

It has been suggested that while this type of assumption is satisfactory 
for isotropic plates, it may not be sufficient for laminated plates and shells 
because of the varying shear moduli through- the- thickness of the laminate. 
Thus, a layerwise constant shear angle theory (LOST) has been used in which, 
as the name implies, a constant shear deformation angle is used for each ply. 

Using the terminology of Bert [33], the majority of composite finite ele- 
ments can be categorized into two distinct groups, namely, (1) a "smeared 
laminate model" (SLM) , and (2) a "discrete layer model" (DLM) . In the SLM 
elements, the laminate is not considered as a series of individual layers but 
a heterogeneous, anisotropic medium. In terms of shear deformation, the con- 
stant shear angle theory (CST) type of elements are in effect smeared models 
since they do not distinguish the shear strain angles on an individual layer 
basis. In the DLM elements, each layer retains its own identity and in some 
cases is treated as a "sub- element" of the total element. In terms of shear 
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deformation, the layerwise constant shear angle theory, (LCST) are discrete 
models since each layer is permitted to have its own shear angle. 

In addition, the above classifications may be further subdivided into 
groups which utilize lower- order and higher- order approaches. The criteria 
for distinguishing between a lower- order versus a higher- order approach de- 
pends on the order of thickness- coordinate (z) terms included in the expansion 
of various displacement components. That is, expansions including only up to 
linear terms in z are described as lower- order approaches, whereas, those in- 
cluding higher- than- linear powers of z are collectively referred to as higher- 
order approaches. It should be noted that application of either approach can 
be made on the entire laminate , as in the SLM elements, or to individual 
layers, as in the DLM element. 

At this point, it may be useful to briefly discuss the higher- order theo- 
ries in general. There have been numerous higher order theories developed, 
Reddy [34], Lo [35] and Murakami [36], to account for shear deformation. In 
these approaches, higher- order expressions containing quadratic and cubic 
terms are used for the displacement or stress field equations. By using such 
expressions, accurate predictions of in- plane displacements, shear strain and 
stress distributions through- the- thickness are obtained. The goal and motiva- 
tion for using these higher- order theories is to be able to obtain an accurate 
prediction of the laminate behavior while still using a 2- dimensional ap- 
proach, as compared to 3- dimensional or elasticity approaches. Thus, the 
higher- order theories attempt to strike a balance between complexity and 
accuracy. 

In the above, two approaches have been used to develop the higher- order 
theories, namely, a displacement- field assumption or a stress- field assump- 
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tion. The most commonly used is the displacement- field assumption, where in 
this case, the new laminate theory is based upon a kinematically admissible 
displacement field. 

The main thrust of the "higher- order theory" models is to primarily deal 
with the complex cross-sectional warping of the laminate. All of the litera- 
ture reviewed, i.e. Reddy. Lo, Murakami, emphasizes how well their particular 
approach can predict the "zig-zag" shape of the in- plane displacement distri- 
butions. In reality, the crucial test appears to be how well the theory pre- 
dicts the through- the- thickness distribution of shear stresses. In fact, 
accurate methods to predict such transverse shear stress variations is cur- 
rently one of the very active research areas in composite modeling; see 
additional related discussions in Sec. 4.2. 

For example, Reddy [34] gives results showing through- the- thickness shear 
stress distributions obtained directly from constitutive equations and equili- 
brium equations. The distributions obtained from the constitutive equations 
are quite poor with stress discontinuities occurring at the layer interfaces. 
The only improvement over the lower- order theory, i.e. CLT, is the parabolic 
shear stress distribution within each layer and the ability to satisfy the 
stress free conditions at the top and bottom surfaces of the laminate. On the 
other hand, if ply equilibrium equations are used to calculate the shear 
stresses, more accurate and correct shear stress distributions are obtained. 
As will be shown later in Sec. 4, the present element can also obtain such a 
prediction through the use of a modified shear strain distribution method, 
which is simpler than using higher- order theory. In the results of Murakami 
[36], the shear stresses are obtained from equilibrium equations, and as ex- 
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pected show good agreement with the exact solution, while Lo [35] makes no 
mention of the shear stress distribution results, only flexural stresses. 

2.1.1 Discrete Layer Model (DLM) Elements 

A group of 8- node isoparametric hybrid- stress elements for thin to thick 
laminated plates have been developed by Spilker [37,38,39]. Together with a 
higher- order approach the hybrid- stress formulation is utilized. The latter 
is claimed to be preferred since it is easier to satisfy interlayer stress 
continuity and laminate surface conditions exactly. The main feature of this 
element is that displacements and non- normal cross- section rotations are inde- 
pendent from layer to layer [37]. Thus, each layer possesses its own degrees 
of freedom (d.o.f.) and remains a discrete layer. In addition, individual 
stress expressions are written for each layer. Stress continuity through- the- 
thickness is then established by relating certain terms of the adjacent z^* 1 
and z'+l th layers. Also, certain displacements of the adjacent layers are 
related to establish displacement continuity through- the- thickness. The DLM 
approach is used in order to account for the severe cross-sectional warping on 
the laminate. This point is emphasized by the fact that within each layer , a 
higher- order displacement assumption is used. Specifically, the in- plane 
displacement, v, is on the order of z 3 and the transverse displacement, w, is 
of z 2 order, where z is the thickness coordinate. This is in contrast to the 
usual applications of the higher- order theory, which is usually applied on the 
laminate level. 

Because each layer is considered to be discrete, forming the element 
stiffness becomes an involved process. Basically, an assembly process is 
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required to form the element, or laminate, stiffness from the individual layer 
stiffnesses. Reference [38] likens the process to the conventional element 
"assembly" operations in finite element codes. That is, stress parameter 
"pointers" and nodal displacement "pointers" are required to add individual 
layer stiffness contributions to the overall element stiffness. Thus, the 
element requires a large amount of computation just to form the element stiff- 
ness. As a matter of fact, Spilker does point out that the number of d.o.f. 
and laminate stress parameters grow rapidly as the number of layers increase, 
and computer core storage may limit the number of layers allowed [37] . 

In reference [38] , in an attempt to eliminate some of the layer depen- 
dency of the element, "kinematic constraints" are introduced. These con- 
straints are such that a number of the individual layer displacements are 
related to a laminate set of d.o.f. located at the mid- plane of the laminate. 
The element d.o.f. are only those associated with the mid- plane reference 
surface making the element’s d.o.f. layer independent. Thus, the element now 
has a total of 40 d.o.f. independent of the number of layers. Also, the 
higher- order displacement assumptions of [37] are simplified to again decrease 
the complexity of the element. Thus, in effect, the detailed in- plane dis- 
placement response is neglected, yet the element still retains the higher- 
order stress expressions. By using such a compromise, the element still gives 
accurate transverse displacement as compared to CLT and good through- the- 
thickness stress distributions as compared to elasticity. 

Finally, in reference [39] Spilker introduces an element which is totally 
layer independent. For moderately thick laminates, 0.01 < h/L < 0.1, shear 
deformation affects structural behavior but the individual layer warping ef- 
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fects are not significant. Thus, individual layer non- normal cross-section 
rotations and individual stress parameters may not be necessary. In this 
element, the total number of d.o.f. and stress parameters are independent of 
the number of layers. First, the nodal d.o.f. are related to the in-plane and 
transverse displacements at the mid- plane reference surface of the laminate, 
and to the non- normal cross-section rotations of the laminate. Second, the 
stress fields in each layer are expressed in terms of laminate stress para- 
meters. All the through- the- thickness integrations can therefore be isolated 
and they are then analytically evaluated. Specifically, this refers to the 
calculation of the laminate material properties, i.e. axial, bending, and 
coupling stiffnesses. In effect, laminate theory is used to calculate these 
quantities, thus effectively eliminating any numerical integration through- 
the- thickness. Since these are the only layer dependent quantities for the 
element, the computation time will increase slowly as the number of layers 
increase. 

The element predicts transverse displacements quite well, as compared to 
exact solutions, for laminate thickness ratios of 0.01 < h/L < 0.25, where h 
is the total thickness of the laminated plate, and L is a characteristic 
"spanwise" length. The stress distributions show good agreement with exact 
solutions for thickness ratios of 0.01 < h/L < 0.10, i.e. up to moderately 
thick laminates. The stress distributions for thick laminates, h/l=0.25, are 
not in good agreement. This is attributed to some of the above mentioned 
simplifications introduced into this element to improve its efficiency. 
Details of improved CPU times can be found in [39], and in general this ele- 
ment proves to be more efficient than the previous elements. 
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Thus, the trend is to move away from the true discrete layer model, which 
is inefficient from a computational standpoint, to a semi- discrete layer 
model. Vhile the DLM elements give very detailed displacement and stress in- 
formation, one needs to question whether or not such detail is necessary. 
Spilker states that such detail may only be required for special cases such as 
non-linear post- first- ply failure, laminate strength measurements, etc. [37]. 
But for typical structural applications such detail is not necessary. 

Chaudhuri and Seide [40] present a triangular 6- node displacement based 
element. It is somewhat similar to Spilker’ s original element [37], in that 
the total d.o.f. of the element is layer dependent. The total laminate ele- 
ment consists of N "layer elements" for a N- layer laminate. Each layer ele- 
ment is triangular in shape and consists of 6- nodes with 5 d.o.f. per node. 
This element is different from Seide ’s previous work in that the element can 
assume a general triangular shape. The reason given for choosing a triangular 
over a quadrilateral shape is that it allows easier modeling of holes and 
cut- outs. 

In terms of the d.o.f., each "layer element" has two independent in- plane 
displacements and two independent rotations. Because each layer has indepen- 
dent rotations, i.e. shear angles, Seide uses the term layerwise constant 
shear angle, 1CST, to describe this type of element. The remaining d.o.f., 
transverse displacement, is not independent from layer to layer. That is, 
each element shares the same transverse displacement with the adjacent nodes 
that have the same in- plane coordinates through- the- thickness. This satisfies 
the assumed condition of transverse inextensibility, i.e. laminate thickness 
does not change, and interlayer continuity. 
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However, Seide’s element differs from that of Spilker’s layer- dependent 
element in one important aspect. That is, because of its underlying displace- 
ment-based formulation, the former element must satisfy the displacement 
continuity requirement, thus allowing the typical numerical integration scheme 
to be used through- the- thickness to form the element stiffness from the indi- 
vidual layer stiffnesses. In contrast, since Spilker used assumed- stress 
fields, which are not continuous through- the- thickness, the complex "pointer" 
method was required. Seide uses Heaviside functions to write a continuous 
function through- the- thickness of the laminate for the in- plane displacements, 
with the resulting expressions being similar to those used in the higher- order 
theory approaches. The expression for transverse displacement is constant, as 
discussed. 

The results produced by this element shows good agreement with the 
elasticity solutions. In particular, the transverse displacement predictions 
for a thick laminate, i.e. L/h = 4, the CLT error is around -787. while the 
LCST and CST give -37. and -5.2 7., respectively. Thus, for thick laminates, 
Seide believes that LCST is required for acceptable results. But it should be 
noted that for moderately thick laminates, L/h = 10 to 20, both the error be- 
tween the LCST and CST theories becomes almost negligible. Thus, since CST 
theory gives adequate results for thin to moderately thick laminates, is LCST 
necessary? 

2.1.2 Smeared Laminate Model (SLM) Elements 

Another classification of composite finite 
termed as smeared laminate model (SML) elements. 


elements are what may be 
Most of the original compos- 
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ite element developments fall under this particular category. The elements of 
Panda and Natarajan [41], Chang and Sawamiphakdi [42], Reddy [43], Rogers and 
Knight [44] , and Haas and Lee [45] , reviewed in this paper, use a "smeared" 
laminate in the sense that the element does not retain the discrete ply struc- 
ture, as do the DLM elements discussed previously. Instead, by either numer- 
ical integration through- the- thickness or the application of laminate theory 
an equivalent laminate or element stiffness is formed. In most cases, the use 
of this smearing technique is sufficient for predicting the global structural 
response. 

The elements of Panda [41] and Chang [42] , can be considered as represen- 
tative examples of early composite finite elements based upon a "degenerated" 
solid shell theory. In particular, the element of Panda is a simplification 
of an element originally developed by Mawemya and Davies [46] . This original 
element considered independent rotations for each layer, thus the total 
degrees of freedom of the element was dependent on the number of layers. 
Panda states that the use of independent layer rotations does not necessarily 
give improved results, but does increase the size of the problem considerably 
[41]. Thus, the element of Panda uses a constant normal rotation through- the- 
thickness of the laminate. Also, details are provided on the method of inte- 
grating through- the- thickness using a "thickness concept" to effectively 
account for the different layer material properties. This change of variable 
scheme is common to most elements which integrate through- the- thickness to 
form the element stiffness. The results provided by Panda show good agreement 
with elasticity solutions. For L/h > 20 errors in deflections and stress are 
on the order of -77. and for L/h > 10, the deflections are within -157. and 
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stresses are within -107.. Panda shows that this particular method of inte- 
grating through- the- thickness provides good accuracy for thin plates and rea- 
sonable accuracy for thick plates [41] . 

A limitation of the previous displacement based "degenerated" elements is 
in dealing with the problem of locking. That is, most shell elements have a 
tendency to overpredict the element stiffness as the thickness becomes small, 
which leads to locking. In reference [45], Haas and Lee present a nine-node 
assumed strain degenerate solid shell element, which was chosen because it has 
been shown that the problem of locking is eliminated. This is achieved by 
segregating the energy expression into components of in-plane, bending, and 
transverse shear strain energy terms. Also, coupling strain energy terms are 
included since coupling effects are present in non- symmetric composite lami- 
nates (for symmetric laminates or isotropic materials, this term will vanish). 
By segregating the in- plane and transverse shear strain terms, the source of 
locking is eliminated. The strain polynomial assumptions were chosen to pro- 
duce a stable element with no spurious kinematic modes [45] . An important 
point is that these strain polynomials are used to describe the behavior of 
the entire laminate through- the- thickness, unlike the stress polynomials of 
Spilker [37]. Because of this, the element degree’s of freedom are not layer 
dependent . 

Results are presented for symmetric laminated plates subjected to a 
sinusoidal load. For square laminated plates with length- to- thickness ratios 
of 10 < L/h < 400, the element shows good agreement with elasticity and thick 
plate theory [45] . One interesting point is that no mention is made of any 
shear correction factors for shear deformation, yet reasonable results are 
obtained. 
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Recall that previously the higher- order theory was shown as applied on an 
individual basis, in particular, the DLM elements of Spilker. On the other 
hand, Reddy [43] uses the higher- order approach on a laminate or element 
level. The element, as described in reference [43], is classified as a SLM 
element where the higher- order displacement fields are used to describe the 
entire laminate through- the- thickness. Specifically, Reddy’s element is based 
upon the previously discussed theory set forth in reference [34]. Again, as 
most proponents of the higher- order theory emphasize, the classical plate 
theory based upon the Kirchhoff-Love assumption is not adequate to describe 
the behavior of moderately thick laminates. By using a higher- order theory, 
the parabolic distribution of the transverse shear stresses is accounted for 
and thereby eliminates the need for the shear correction factors, [43]. 

Results are presented covering deflections, stresses, buckling loads, and 
vibration frequencies. If attention is focused primarily on deflection and 
stress results, the higher- order element does show excellent agreement with 
exact elasticity solutions for the cases presented here. Reddy again empha- 
sizes the discrepancy between the CLT and HSDT, with the difference becoming 
more significant as the ratio of longitudinal to transverse moduli increases. 

But as is common with the higher- order approach, Reddy does point out the 
associated drawbacks. First, higher- order stress and moment resultants arise 
from the formulation which are difficult to understand physically. This same 
problem is evident in the elements of Spilker. And secondly, the theory con- 
tains derivatives of the transverse deflections. Thus, the element requires 
functions with C 1 continuity [43] . 
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In reference [47] , Lakshminarayana and Murthy present a shear flexible 

i 

finite element which is based upon a generalized theory of Yang, Norris, and 
Stavsky, YNS. It is "shear flexible" in the sense that shear deformation is 
accounted for. Using this theory, the coupled bending and stretching response 
of the laminated plate is expressed in terms of five unknowns. That is, three 
displacements; u,v,w, and two rotations; e x ,© y , of the normal [47]. The 
element itself is triangular in form and has three nodes located at the ver- 
tices. Complete cubic polynomials are used to approximate the displacements 
and rotations within the element and were chosen because the authors believe 
that higher- order elements have advantages over lower- order elements when 
dealing with laminated composite plates and shells [47] . There are 15 
degrees- of- freedom per node; three displacements, two rotations, and their 
derivatives. The authors believe that by using the derivatives of displace- 
ments and rotations at the nodes is an advantage because it allows evaluation 
of the stresses directly at the nodes instead of at Gauss points which is 
common in most other elements. But this also makes the element very computa- 
tionally expensive since for three nodes with 15 degrees- of- freedom per node 
gives a total of 45 degrees- of- freedom for the element. Thus, one must ques- 
tion whether the luxury of having the stresses directly at the nodes is really 
advantageous. Again, as with previous SML elements, the element stiffness is 
partially based upon laminate theory to reduce the multi-layer composite into 
a smeared single layer plate. 

Lakshminarayana and Murthy present a large number of cases with a variety 
of geometry, loading, support conditions, and plate thicknesses to access the 
performance of the element and the importance of shear deformation effects. 
In each case, the performance is measured by comparing to either analytical or 
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other finite element solutions obtained from a reference by Noor and Mathers. 
The element appears to agree very well with the other solutions, both finite 
element and exact, in all cases presented. Because of the large amount of 
detailed results and the varied cases, this reference will serve well as a 
source of benchmark problems for the present composite finite element 
development. 

Rogers and Knight [44] present a higher- order element with a specific 
application to filament- wound composite pressure vessels. It is pointed out 
in the paper, that the main motivation for the development of the element is 
to provide an efficient and simple means of modeling the complex composite 
structure found in the pressure vessels. Also, because of the low shear 
stiffness of the material, significant shear effects between the layers 
develop, thus requiring a higher- order displacement field. The proposed ele- 
ment is similar to the higher- order elements previously discussed, but a dis- 
tinguishing feature presented is that the order of the through- the- thickness 
polynomial is user selectable. The element subroutine internally generates 
the appropriate shape functions using Lagrange interpolation polynomials. 
This is unlike most elements in which the shape functions are permanently 
coded into the element subroutines. This eliminates the need to manually 
derive the shape functions and their derivatives. Whether or not this feature 
is really an advantage is not clear. By using such a higher-order element, 
the typical multi-layer composite can now be effectively modeled using one 
element through- the- thickness. As proof of the need of a higher- order ele- 
ment, data is presented which shows the amount of computation time saved by 
using one higher- order element through- the- thickness, as compared to the use 
of individual elements for each ply. 
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The element stiffness is generated by performing numerical integration 
through- the- thickness on a layer-by-layer basis. Unlike most composite ele- 
ments reviewed that use this technique, this element uses a zeroth- order rule 
(rectangular integration rule) with the weighting factors being proportional 
to the layer thickness in which the integration point resides. This allows 
the effects of each individual layer’s material properties to be effectively 
modeled. The authors do point out though, that no relationship exists between 
the polynomial order, n, and the number of integration points required to 
integrate the polynomial exactly using this zeroth- order integration (unlike 
standard Gauss quadrature) . This creates a variable in the accuracy of the 
finite element solution which can be only eliminated by experience with this 
particular element [44] . 

All of the papers reviewed thus far, have been limited to linear problems 
with regard to deflections, vibrations, and buckling loads. One of the ob- 
vious reasons for such limited work in nonlinear composite analysis is the 
added complexity of the problem. Indeed, considering composite literature in 
general, closed- form or "exact" solutions that do exist are further limited to 
a very narrow class of problems; i.e., for the most part, to symmetrically 
layered cross- ply and angle- ply laminates. This is due to the presence of 
coupling terms, which are peculiar to composite materials, that often make any 
analytical solution for unsymmetric laminates intractable. This accounts for 
the repeated use of the 0/90 class and other symmetric laminates, which do not 
contain these coupling terms, as reference solutions in most finite element 
literature; e.g. Pagano [48] and Vhitney [49]. 
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A comprehensive review of research on both linear and nonlinear analysis 
of composites was given in [50]. In addition, a penalty- type formulation to 
overcome shear locking for thin plates together with the composite theory by 
Yang, Norris, and Stavsky (YNS) , and various classical first- approximation 
shell theories were utilized in [50] to develop a displacement- based element 
for nonlinear analysis. Similar applications of these types of elements can 
also be found in [52,53]. 

3. Present Element Formulation 

The element to be presented here is the extension of earlier work by 
Saleeb et.al. [27] to the modeling of multi-layered composite materials. It 
is the objective of this investigation to see how well the element will per- 
form with the minimum amount of modifications to its basic formulation. If 
the element yields satisfactory results, as such, a good foundation will have 
been established in which to build upon for future work. It is also desirable 
to keep the element as "simple” as possible thereby minimizing the computa- 
tional expense, i.e. an element which has the total d.o.f. independent of the 
number of layers. 

This element may be classed as a hybrid/mixed or assumed- strain shell 
element. That is, in the underlying variational principle, both the displace- 
ment and strain fields are approximated independently. Since the ultimate 
goal of the research is the application of the elements for high- temperature 
metal- matrix composite structural analysis, i.e. non-linear material analysis, 
it appears that utilizing a strain assumption, as opposed to a stress assump- 
tion, is most desirable. That is, in the case of material non-linearity it is 
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still entirely valid to assume a linear variation of bending strains through- 
the- thickness of the "thin" laminate, while such an assumption for associated 
stresses is certainly not valid. Also, choosing a strain field a priori is 
much easier to justify than it would be for a particular stress field assump- 
tion. 

Vhat is different in the present approach from some of the previously 
mentioned elements is that the assumed strain field used represents the entire 
laminate. Thus this element can be classified as a smeared laminate model, 
SML, type, utilizing a first- order or constant shear angle (CST) theory. Its 
general formulation is derived from a "modified" form of the Hellinger- 
Reissner variational principle, as described in the sequel. 


3.1 Mixed Variational Principal 

The "modified" Hellinger-Reissner principle, including the kinetic energy 
term for the dynamics problem, is expressed as 
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where 


e = independently assumed strains 
e = B^q = strains derived from assumed displacements 
D = material stiffness matrix 
p = material density 
u = assumed element displacements 
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u = element velocities 
T = boundary tractions 
f = boundary force 
V = volume 
r = time 

All quantities with an overbar, ("), are prescribed quantities, and and S u 
are portions of the element boundary over which tractions and displacements 
are prescribed. By taking the first variation of the functional with respect 
to the strains e, and displacements u, and invoking the stationarity 
conditions = 0, 


s ’ti - 4Pvj 
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dr (2) 


the equilibrium equations, strain- displacement relations, i.e., 
e = e = B u in V 
B^<r + f = pii in V 

together with the associated surface- tractions and displacement boundary 

conditions on and S u , respectively, are recovered. In the above, B is a 

T 

differential strain- displacement operator, and B is the associated "matrix- 
divergence" operator. 


3.2 Element Stiffness and Mass Matrices 

The derivation of the stiffness and mass matrices will be obtained from 
the previously described functional r^, now referring to lamina- coordinate 
systems. Justification for choosing such a reference frame is given in detail 
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in our previous work [27]. Thus, the functional can be expressed as, 
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where dV = dxdydz = | J | drdsdt , in which r,s,t are the natural (isoparametric) 
coordinates and J is the Jacobian- of- coordinate- transform matrix. For the 
finite element discretization, the displacements, u, are interpolated from the 
nodal degrees- of- freedom, 
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where the are the two-dimensional shape functions in terms of local (iso- 
parametric) coordinates, r and s, for node k. For simplicity in terms of the 
derivations, the above equation can be written symbolically as, 


( 5 ) 

where N are "modified" shape functions obtained by combining terms in Eq. (4). 
In the above, e^®\ with i = 1,2,3, are the components of the fiber triad 
associated with node (k), which are defined as in [27]. 

The independently- assumed strains, e, are approximated in terms of the 
strain parameters, /?, as follows 

£ = P § (6) 

where P is the strain interpolation matrix. Again for brevity, we refer to 


u = Nq 
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[27] for details of the form for the interpolation matrix P. The strains, 
derived from the assumed displacements are expressed as 


r\ 

e 



( 7 ) 


l 

where B is now a local strain- displacement matrix. 

Now, substituting Eqs. (5), (6) and (7) into Eq. (3), the functional 
becomes 
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Equation (8) can be simplified and written as 
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where the matrices H, G, H and Q are defined as, 
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In order to eliminate the strains on the element level, stationarity is 
invoked on Eq, (8) in terms of (3 variables, thus yielding 


g = H' hq 


(14) 


Now, substituting Eq. (14) back into the functional and taking the 
variation with respect to nodal displacement parameters, we finally get the 
element stiffness equations as follows, 

(G T H _1 Gq + Mq - Q)£q T = 0 
or 

Mq + Kq = q (15) 


where the stiffness matrix K is defined as 

! = g 1 !' 1 ? ( 16 ) 

and the mass matrix M is defined previously in Eq. (12). 

3.3 Through- the- Thickness Integration for Composites 

According to Eqs. (10) and (11), the component matrices for the present 
"smeared" element stiffness (see Eq. 16) are obtained by simply integrating 
through- the- thickness and taking into account the differences in material 
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properties of each layer. That is, since the element now represents a multi- 
layered composite, the material matrix D is no longer constant through- the- 
thickness. Thus, individual layer material matrices, D^, must be formed, 
transformed, and then "combined" during integration to form the smeared 
laminate stiffness. 

Specifically, each layer is assumed to be in a state of plane stress 
(< 7-33 = 0 ) , i.e. the stresses in local material fiber or symmetry- axes coordi- 
nates are ordered as 

{ ff )l ~ { *11**22’ r 12 ,r 23 ,r 13 


and the corresponding layer constitutive relation is, 
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or symbolically, 


{*} i ~ 


(18a) 


(18b) 


in which the matrix elements are defined as, 
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where and E ^2 are the longitudinal and transverse elastic moduli, re- 
spectively, the G’s are the transverse and in- plane shear moduli and i/’s are 
Poisson’s ratios, with all quantities being defined in material fiber coor- 
dinates. Now, the layer matrix needs to be transformed from material fiber 
direction to the direction of the global axes. This is accomplished by using 
the standard coordinate- transformation relation: 


m* = [»]< m* 


where [D] ^ is the transformed layer material matrix and [T] is an appropriate 
transformation matrix. This latter matrix will generally vary along the 
lamina surface (i.e., t=0 surface) from one integration point to another (with 
different r,s- coordinates) , depending on the composite’s material- fiber 
layout, in a generally- distorted (non- rectangular) and curved shell element. 
However, for the special case of flat plate problems with "regular" meshes, as 
in the numerical simulations reported later in Sec. 4, a simplified unique 
transformation matrix can be utilized; i.e., 


[T] 


£ ~ 


cos 2 © sin 2 © 2sin©cos© 0 0 

sin 2 © cos 2 © -2sin©cos© 0 0 

-sinecose sinecos© cos 2 ©- sin 2 © 0 0 

000 cos© sin© 

000 -sin© cos© 


( 20 ) 


where the quantity © is the angle the material fiber makes with respect to the 
x-axis, as shown in Fig. 2. 

Thus, the stress- strain relation in "global" coordinates for layer £ is 
finally written as 


35 



( 21 ) 


f * 


Pll 

D12 

Di 3 

0 

0 


Cx ' 

'y 


p2 1 

I>22 

D23 

0 

0 


€ y 

Txy 

• = 

D3I 

®32 

D33 

0 

0 


7 xy 

Tyz 


0 

0 

0 

&44 

D45 


7 yz 

. T xz. 

i 

_ 0 

0 

0 

D54 

D55 

l 

7 x Zj 


Finally, in order to facilitate the integrations indicated in the above H 
and G matrices, we utilize the following simple "change- of- variable" procedure 

[41]: 
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which yields, upon differentiating with respect to t^, 
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(See Fig. 3 for definitions of the various quantities used in the above 
expressions). Thus, the matrices H and G now have the forms, 
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which are evaluated using standard numerical integration; i.e., two Gauss 
quadrature points per layer, and the indicated summation is over the total 
number of layers, n. 
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4. Numerical Studies 


In this section, we present a series of numerical simulations for select- 
ed problems taken from various literature sources, to demonstrate the capa- 
bilities of the element in representing the static and dynamic behaviors of 
composite laminates. These problems involve varied boundary conditions, 
aspect ratios, loading conditions and laminate configurations. For several of 
these problems, comparisons are made with available "exact" analytical solu- 
tions as well as other independent finite element solutions. The quantities 
of interest are predictions of laminate global response, i.e. deflections, and 
natural frequencies with corresponding vibration modes, and laminate local 
responses, namely, bending, in- plane, and transverse shear stresses. For 
transverse shear stresses in particular, a series of alternative methods are 
investigated in an attempt to "improve" its accuracy. 

In the following results, the elements of Lakshminarayana and Murthy will 
be conveniently referred to as TRIPLT [47] and those of Spilker as MQH3R [39] 
and V2R [28] . These particular elements are chosen because they appear to 
provide accurate results for comparison. 

4.1 Static Test Problems: Deflections and Stresses 

4.1.1 A Cantilever Beam with In- Plane Point Load 

The following problem, taken from Lakshminarayana and Murthy [47] , is 
used to access the element’s prediction of the in- plane response of a lamina- 
ted composite beam; i.e., a slender cantilever beam subjected to an in- plane 
point load of magnitude P = 1.0 lb (Fig. 4). The material used has = 30 x 
10® psi, E 22 = 0.75 x 10® psi, = 0.45 x 10® psi, v ^ = 0.25. A variety of 
laminate configurations are used, i.e. 1- , 4-, and 8- layer laminates. Such a 
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variety of multi-layer laminates will provide for a validation test of the new 
element subroutines for layer- transformation and integration through- the- 
thickness schemes. The "exact” solution is based upon that of Lekhnitskii, 
and the TRIPLT element is used for comparison. 

As shown in the results summarized in Table 1, the present element com- 
pares quite well with the exact solution in all cases analyzed. In addition, 
although its "overall" accuracy is comparable to that of TRIPLT, the present 
element’s results is evidently superior in the e = 30 single- layer case. 

4.1.2 A Single Layer Clamped Rectangular Plate Under Uniform Pressure 

The problem analyzed here is a single layer rectangular plate subjected 
to a uniformly distributed pressure load of magnitude q Q . This problem was 
also taken from Lakshminarayana and Murthy, and even though no "exact" solu- 
tion is available, it is felt that the TRIPLT element has proven to provide 
accurate results. This problem is a difficult test for the element due to the 
combination of plate aspect ratio, boundary conditions, and high material 
anisotropy. The geometry of the plate has an aspect ratio of 2, and has to- 
tally clamped boundary conditions, Fig. 5. The material used has the follow- 
ing properties, E^ = 30 x 10®, E 22 = 0.75 x 10®, = 0.45 x 10®, Ggg = 
0.375 x 10® u n = 0.25. 

Because of the demanding nature of the problem, a mesh convergence study 
was performed to assess the elements behavior as a function of mesh size. Due 
to a lack of material symmetry, the entire plate was modeled. The results for 
the TRIPLT element were obtained using an 8 x 8 mesh. Recall the TRIPLT ele- 
ment uses 3 nodes with 15 d.o.f. (degrees of freedom) per node, thus giving a 
total of 1215 d.o.f. for the 8x8 mesh. The present element’s 10 x 10 mesh 
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contains 605 d.o.f. As shown in the convergence results in Table 2, for most 
material fiber angles, convergence is attained at the 8x8 mesh with the 10 x 
10 mesh giving the best results. Thus, the present element is able to model 
the problem, accurately using less d.o.f. 

Table 3 summarizes the above results for normalized center deflections in 
a 10 x 10 mesh for the present element and an 8 x 8 mesh for the TRIPLT 
element. 

4.1.3 Two- Layer Clamped and Simply Supported Square Plates Under Uniform 
Pressure 

A square, 2- layer plate subjected to a uniformly distributed pressure p Q 
is considered here. The plate was analyzed using both simply supported and 
clamped boundary conditions. Note that the simply supported boundary condi- 
tions have, in addition to the transverse displacement, the in- plane displace- 
ment normal to the plate edge restrained. The clamped boundary conditions 
have all d.o.f. restrained along the plate edge. The material properties for 
the problem are; = 40 x 10® psi, E 22 = 1 x 10® psi, = 0.5 x 10® psi, 
v^2 ~ 0-25. 

The results are compared to exact solutions as well as other finite ele- 
ment results using the TRIPLT element, and the MQH3T and V2R thin plate 
elements of Spilker [38,39]. Both of these latter elements contain 8 nodes 
with 5 d.o.f. per node. A 6 x 6 mesh of MQH3T elements is used for the clamp- 
ed plate (665 d.o.f.) and a 4 x 4 mesh of V2R elements is used in the simply 
supported plate (325 d.o.f.). The TRIPLT element utilized a 6 x 6 mesh (735 
d.o.f.). The present element used a 10 x 10 mesh for both cases (605 d.o.f.). 
Again, due to a lack of material symmetry, the entire plate was modeled, and 
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the results are summarized in Tables 4 and 5. Clearly, the present element 

agrees well in both cases with the exact solutions. 

In addition to displacements, bending moments M and H were calculated 

xx yy 

for the simply supported plate, and are given in Tables 6 and 7. 

4.1.4 Cylindrical Bending of Symmetric Laminates Under Sinusoidal Pressure 
In the present example, the cylindrical bending of a symmetric laminate 
is considered. The problem is a simply supported plate subjected to a sinu- 
soidal pressure load of maximum intensity q Q . The material properties are; 

= 25 x 10® psi, Eg 2 = 1 x 10® psi, = 0.5 x 10® psi, G 22 = 0.2 x 10® 
psi, z/j 2 = 0.25. 

Exact solutions for this specific problem have been obtained by Pagano 
[56,57], providing detailed stress distributions, in addition to deflections. 
Thus, this example will provide a means of determining the present elements 
prediction of laminate local behavior; in particular, transverse shear stress- 
es, which will be discussed in the next section in some detail. 

Figures 6 and 7 show deflections, w, as predicted by the present element 
and compared to the exact solutions, where ' 
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Figures 8 and 9 show the predicted through- the- thickness, (T-T-T) distri- 
bution for in-plane stress, , which matches that of Classical Plate Theory 
(CPT) quite well, as expected. 
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4.2 Treatment of Transverse Shear Stress 

As mentioned previously, one of the desirable features of any composite 
element is its ability to accurately predict the "detailed 11 stress distribu- 
tion within the composite. However, in this regard, one of the obvious limi- 
tations of the element presented here is the assumption of constant transverse 
shear strain through- the- thickness (see Fig. 10 with the same problem as des- 
cribed in Section 4.1.4 being considered here). The resulting transverse 
shear stress, if directly obtained from the constitutive relations, are dis- 
continuous through- the- thickness, and as Fig. 11 shows, will be grossly in 
error. Thus, some consistent means needs to be established to allow more 
accurate transverse shear stress distributions through- the- thickness. Again, 
it was still desired to avoid any modifications to the element’s basic formu- 
lation itself. Thus some form of "post- processing" of the transverse shear 
stresses was attempted while still keeping the "computationally- efficient" 
global scheme of utilizing constant through- the- thickness and layer- number 
independent distributions for transverse shear strains . To this end, several 
approaches are described next. 


4.2.1 Approach 1: Modified Shear Stress Distribution 

The first approach involves the simple redistribution of the transverse 
shear strain into a parabolic form through- the- thickness. This type of 
assumption is believed to be more reasonable in that the conditions of zero 
shear strain at the top and bottom surfaces are automatically satisfied. The 
value of the constant shear strain obtained directly from the finite element 
results 7a ve » is redistributed parabolically, i.e. 


'max 


It 


ave 


7(z) = 7, 


max 



[„•) 2 

1 + 4 

h] . 


(25) 
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The corresponding transverse shear stress is then computed using "layer" 
constitutive equations, 


{r} k = [D] k { 7 } k (26) 

where the subscript k corresponds to layer k. 

As shown in Fig. 12, using this method does produce more accurate trans- 
verse shear stress distributions, yet the stresses remain discontinuous at the 
layer interfaces. Note, however, that the magnitudes of the shear stress are 
closer to the exact values. It may also be noted that this type of distribu- 
tion is similar to that of a higher- order theory presented by Reddy [34]. 

4.2.2 Approach 2: Equilibrium- Based Method 

In the alternative method investigated here, an attempt is made to elimi- 
nate the transverse shear stress discontinuities noted previously at the layer 
interfaces. In particular, the methodology outlined by Chanduri and Seide 
[58,59] was chosen. In this approach, "layervise" equilibrium and stress con- 
tinuity conditions are imposed. One problem arises, however, in choosing the 
values of average layer shear stress, r k - Note that the exact methodology as 
given in [58,59] cannot be applied here to achieve the same good results due 
to the rather significant differences between the elements utilized. As dis- 
cussed previously, the element of Seide’s, for which this method was devel- 
oped, is a layer- dependent element of the displacement type, i.e. the total 
degrees of freedom are layer- dependent quantities. In contrast, the mixed 
element presented here, has in effect, only one "layer" to model the entire 
laminate thickness. 
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Initially, the average shear stress values used were those obtained di- 
rectly from the finite element analysis (refer to Fig. 11). The transverse 
shear stress distribution as compared to the exact 3-D solution is shown in 
Fig. 13. The shear stresses in the middle 90- layer are grossly underestimated 
and the stresses in the outer two layers (both 0) are overestimated. Obvious- 
ly, no improvement has been attained. Thus, two additional cases using dif- 
ferent values for the average layer shear stress were investigated. In one 
case, the average layer transverse shear stresses were calculated by taking 
the average distribution as predicted by the previous Approach 1; see Fig. 12. 
Using these values for r^, produced the distribution shown in Fig. 14, where 
we note that the shear stress values in the middle layer are still underesti- 
mated. In the second case, the value of was taken to be the value of the 
average shear stresses at the interface of any two adjacent layers, and the 
resulting distribution is shown in Fig. 15. The persistent problem of the 
poor stress prediction in the middle layer remains. But it should be noted 
that at the layer interfaces the value of the shear stresses match those of 
the exact solution rather well. 

4.2.3 Approach 3: Simplified Strength- of- Materials Concept 

Finally, the third approach investigated is based upon a basic strength 
of materials concept, specifically, the transformation of the various layers 
to "equivalent" cross sectional areas depending on their elastic moduli. This 
is a method used frequently in the design of heterogeneous materials (e.g., 
reinforced concrete). The layer cross-sectional areas are modified by using 
the following "modular" ratio, n, 
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(27) 


n 



where E is the selected "reference" elastic modulus and E, is the elastic 
0 & 

modulus for layer k. Thus the equivalent cross-sectional area is simply, 


/ 




(28) 


Once the cross section of the laminate is transformed into an equivalent 
section, the transverse shear stress distribution is then calculated by the 
following formula, 


r ( z ) - if ( 29 ) 

where V is the total shear force acting on the cross section; Q the statical 
moment of cross- section about the neutral axis; I the moment of inertia of the 
transformed cross- section; and t the width of the layer at a point z. 

The calculation of the total shear force V, per unit length, is obtained 
from the transverse shear stress distribution obtained directly from the 
finite element solution (Fig. 11), i.e., 


V = S r.h. (30) 

k=l K K 


The results of this method are shown in Fig. 16, and clearly illustrate 
the excellent stress distribution obtained. The quality is far better than 
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either of Approaches 1 or 2. Of course, the major limitation to this approach 
is the fact that the expression for transverse shear is based upon a linear 
distribution of stress, thus linear material behavior is necessary. Since the 
immediate goal is to extend the current element formulation to non-linear 
composite material behavior, this type of approach will need to be modified. 


4.2.4 Additional Remarks on Alternative Approaches for Transverse Shear 
Stress Computations 

Approaches 1 and 2 above produced rather "mixed" results. Approach 1 is 
still attractive from the standpoint that the assumptions made for the current 
linear analysis are still applicable to non-linear analysis. But, as mention- 
ed, the drawback is the stress discontinuities at the layer interfaces. Con- 
cerning Approach 2, a point of difficulty is in the determination of the 
values for r^. A somewhat "ad hoc" procedure was used and such a methodology 
is not desirable. 

With regards to Approach 3, it is certainly appealing in view of its 
extreme simplicity. Some initial thoughts have already been given for the 
extension to non-linear analysis, and it appears that some workable extension 
is possible. In this, the restrictive condition of linear through- the- thick- 
ness distributions of in- plane normal stresses will be replaced by a modified 
"piecewise- linear" assumption for each layer, and the "local" form of equili- 
brium conditions are imposed to derive the applicable formula for transverse 
shear stress distribution. Whether or not this approach will work as well and 
provide as accurate results in the non-linear regime as it did for the linear 
case, remains to be seen. Some additional numerical studies in the linear, as 
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well as the non-linear cases will need to be conducted. And this is a part of 
our planned future research work. 

The difficulties associated with accurate distributions of transverse 
shear stresses serve to show the inherent limitations of the "smeared" com- 
posite element approach. The accurate distributions obtained by other compo- 
site elements are at the expense of using rather computationally- expensive 
formulations. Very detailed distributions are simply not possible using the 
type of elements as the one presented here. But again, one must question how 
crucial is it that the element be capable of detailed "local" quantities in 
addition to accurate "global" quantities, such as displacements, vibration 
frequencies, etc. It is our view that such detailed "local" responses will 
always necessitate very refined, and also complicated, formulations that can 
best be used only for critical parts of the laminated structure determined on 
the basis of a less- expensive "overall" analysis of the "global" response. 

4.3 Test Problems for Dynamic Analysis 

As described in [27] , a feature of the element presented here is that the 
element is defined by five nodes, four of which are defined externally and the 
fifth node being generated internally. The fifth node’s degrees- of- freedom 
are then condensed out of the element stiffness matrix before assembly into 
the global stiffness matrix. It is this feature which presents some difficul- 
ty when dealing with the mass matrix for dynamic analysis. This point can be 
made more clear by looking at the eigenvalue problem being solved. That is, 
the associated eigenvalue problem can be expressed as, 


(K - AM) {d} = 0 


(31) 
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where K and M have dimensions 25 x 25. In matrix form, 


'r K u 
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dc 
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( 32 ) 


For the case of lumped mass, which is being considered here, Mi 2 and M 2 i 
are zero. Thus, the eigenvalue problem is reduced to, 
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dm 


= 0 


(33a) 


Using simple algebra and solving for middle node degree- of- freedom, d„,, 
we have, 

(33b) 


dm = - (K 22 - AM 22 )-% 2 T d c 


and the corner node degrees- of- freedom, d c , are 


(K t 1 - AM 11 ) d c + K 1 2 d m = 0 (33c) 

Substituting the expression for dm into the above, we have, 

(Kn - AM n )d c + Ki 2 (- (K 22 - AM 22 )* 1 K 12 T d c ) = 0 (33d) 

As is evident in the above expression, the equation is coupled containing both 
known, stiffness and mass, terms, and unknown, A, quantities. Particularly 
the term (K 22 - AM^)* 1 presents computational difficulties as discussed by 

Kidder [60] . Ordinarily, some iterative process would be required to solve 
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the above. Most of the literature dealing with condensation in general, 
[ 61 - 63 ], has been applied to the global mass and stiffness matrices. That is, 
the focus has been in reducing extra structural degrees- of- freedom through the 
use of condensation techniques. Thus, iterative methods are acceptable on the 
global mass and stiffness matrices. But in the present case, the condensation 
is required on the element level. Thus, the above methods would have to be 
applied to each element which would be computationally expensive. 

With the above discussion in mind, a suitable alternative to condensation 
is desired. One such alternative is by neglecting the middle node entirely, 
thereby eliminating the need for any form of condensation. Specifically, the 
lumped mass matrix M, as defined in Eq. 12 of Sec. 3.2, is constructed using 
shape functions, N which are based upon four boundary nodes only. By using 
the four- node shape functions, the mass will be "lumped" at the corner nodes 
only. It is believed that this approach is justifiable for the following 
reason. The middle node’s additional displacement is very small relative to 
an "average" displacement interpolated from the four corner nodes. Thus by 
using a sufficiently fine mesh, the contribution of the middle node becomes of 
less importance. Therefore, the current assumption is that neglecting the 
middle node may not lead to significant errors in frequencies. In the remain- 
der of this section, a series of numerical problems are investigated to deter- 
mine the validity of the approach outlined above. Since this is the first 
implementation of the dynamic analysis capability for the element, both iso- 
tropic and anisotropic applications are required to test the new subroutines 
for dynamic analysis. 
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4 . 3.1 Isotropic Material Applications 

4 . 3 . 1.1 Isotropic Cantilevered Plate 

This first example of an isotropic cantilevered plate is chosen to pro- 
vide a method of determining how significant neglecting the middle node is 
over range of frequencies. This particular problem, found in reference [65], 
provides both numerical and test results for the first twelve frequencies of 
vibration. In the analysis performed here, mesh sizes that range from 4x4 to 
10 x 10 , for the entire plate, are used to study the effect of mesh refinement 
on the convergence of the first six frequencies. 

As Table 8 shows, reasonable results are obtained for the first two 
frequencies for all mesh sizes considered, (i.e. the error is within 87. for 
6 / 2 ). When considering the higher modes, as shown, additional mesh refinement 
is necessary in order to produce reasonable results, (i.e. error within 77. for 
u 6 ). Thus it appears that for the first six frequencies, the element gives 
accurate results. 

4.3.2 Composite Material Applications 

In this section, a series of composite structures are examined with the 
specific examples chosen to study the effects of fiber orientation, number of 
layers, aspect ratio, and boundary conditions on the frequency response. In 
most cases, due to the lack of "exact" solutions, comparisons are made to 
other finite element results. Specifically, comparisons are made to the 
displacement- based elements of Oden [24], Owen [64], and Reddy [43,51]. In 
all of the following examples, one quarter of the plate is modeled using a 6x6 
mesh. 


49 



4.3.2. 1 Anisotropic Cantilever Beam 

A simple cantilevered beam which has a single layer with varying fiber 
orientations is considered here. The reference "solution" is that provided by 
Oden et.al. [24]. The material properties used are those which are said to 
represent cubic synonymy that simulates the single crystal structure of nickel 
alloy found in turbine blades [24] . Thus in this case, the fiber angle is 
actually the preferred direction of the single crystal material and is allowed 
to vary from 0 to 90, Fig. 17. Specifically, the material properties are 
Ei i -E 2 2 =1 - 9716 x 10 6 psi, Gi 2 =5.4758 x 10 6 psi, i/i 2 =0.2875, and p= 0.25. In the 
analysis of Oden, both linear and quadratic two-dimensional elements are used. 
Namely, the analysis uses a mesh of 2x40 linear elements (I) or a mesh of 1x10 
quadratic elements (II). 

Tables 9 and 10 show the results of our analysis as compared to that of 
Oden et.al. for the first two frequencies. It should be noted that the ele- 
ment mesh used for this analysis was constructed so that the vibration dis- 
placement is out- of- plane of the element as compared to Oden where the dis- 
placement is in- plane. Thus, a "simpler" mesh of 1x20 elements is used. As a 
matter of fact, a mesh of only 1x10 elements gave satisfactory results, but a 
slightly finer mesh is used so that plots of the corresponding mode shapes 
will be more accurate. Fig. 18 shows the mode shapes corresponding to the 
first two frequencies of vibration for a crystal direction e = 0°. The mode 
shapes compare well to those given in [24] . 

4. 3. 2. 2 Four Layer Symmetric and Anti- Symmetric Square Plates 

This example is that taken from Reddy [51] . Two different laminate 
configurations are considered: 1) a [45/- 45/- 45/45] symmetric laminate, 2) a 
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[45/- 45/45/- 45] anti- symmetric laminate. Both simply supported and clamped 
boundary conditions are used as described in [51] . Two different materials 
are considered: Material I: Eu=25 x 10 8 psi, E 22 =l x 10 8 psi, Gi 2 = 0.5 x 
10 6 psi, G 23 = 0.2 x 10 8 psi, i'i 2 =0 . 25 . Material II: En=25 x 10 8 psi, E 22 =l x 
10 6 ps i , G 1 2 =0 - 5 x 10 8 psi, G 2 3 =0 . 2 x 10 8 psi, 1 ^ 12 = 0 . 25, with p=l for both 
material types. 

Fig. 19 shows the comparison of the present element with that of Reddy’s 
and with closed- form solutions where noted. As is shown, the present element 
shows good agreement in all cases. It even shows better accuracy than Reddy, 
when compared to the closed- form solution for the simply supported material II 
case. 


4. 3. 2 . 3 Simply Supported Square Laminated Plates with Varying Number of 
Layers 

The present example is chosen to demonstrate the effect of thickness and 
the number of layers on the fundamental frequency. A number of cases were 
investigated, namely; two, three, and eight- layer anti- symmetric laminates. 
In all cases, layers of the same fiber orientation have the same thickness 
while the total laminate thickness remains constant. The plates have simply 
supported boundary conditions, that is, in- plane displacements are fixed in 
the tangential direction and symmetry boundary conditions enforced since only 
one quarter of the plate is modeled. Again, a 6 x 6 mesh is used in the present 
analysis. The material properties are the same as those described as Material 
II in section 4.3.2. 1. Fig. 20 shows the results for the three- layer [0/90/0] 
laminate, where comparison is made to both the finite element solutions of 
Owen [64] and classical plate theory, CPT. 
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As shown, for large laminate thickness, the difference between CPT and 
the finite element solutions is significant, and therefore, CPT is inaccurate 
for thick to moderately thick laminates. Again, the pronounced difference 
between the CPT and finite element solutions is noted. In both figures, the 
present element agrees well with the results of Owen. Unfortunately for 
comparison purposes, Owen makes no mention of the mesh size used in the 
analysis. But relatively speaking, these simple problems serve to demonstrate 
the present element’s ability to produce accurate results. 

4.3.3 Additional Remarks 

Considering the results presented above, it appears that the current 
assumption for forming the element mass may be acceptable. Of course, no firm 
conclusions can be made at this time because of the limited, and somewhat 
simplistic numerical studies performed thus far, that is, only simple beam and 
flat plate structures are considered. The problems of the above type have 
predominantly translatory mass effects. The rotary inertia terms have yet to 
be examined. Thus, additional studies, in which curved shell geometries are 
modeled, need to be investigated. 

5. Final Conclusions 

In this part of the report, the formulation and numerical performance of 
the mixed element as extended to laminated composite applications is presented 
with both linear static and dynamic analysis capabilities demonstrated. Even 
though some aspects of the element require further investigation, as outlined 
in Sec. 4.2.4 and 4.3.3, a certain degree of confidence is developed consider- 
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ing the encouraging numerical results obtained in this report. The immediate 
goals of the current research are as follows. 

The extension to shell problems is the next immediate objective. At the 
present, the only extensions that are foreseen in order to make the element 
applicable to shell problems is the generalization of the transformation 
matrix. That is, a full three-dimensional transformation matrix will need to 
be used. Such a transformation matrix already exists in in the computer code 
utilized (NFAP) , thus making the extension straight-forward. Additional test 
problems of shell type applications will be made. 

Some thought has been given to possible refinements for the shear correc- 
tion factor /c. As is consistent with the Reissner-Mindlin type of elements, a 
shear correction factor, k, is chosen to account for shear deformation. It 
should be noted that the value of the shear correction factor currently used 
is 5/6 and is typically used for isotropic cases. It may be more appropriate 
to select a value for k which depends on the laminate configuration being 
analyzed. A methodology for this as presented by Vhitney [54] and Noor [55] 
has also recently proposed a method of an a posteriori correction of the k 
values. Vhether or not such modifications are necessary remains to be seen. 

And finally, the extension to non-linear material behavior will be the 
primary focus of the research. Modifications to the element formulation to 
account for this type of analysis will be made. Part IV of this report will 
outline this aspect in more detail. In general, it is believed a firm basis 
is established in which future research may build upon. 
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PART Ills COMPUTER IMPLEMENTATION AND NUMERICAL METHODS FOR VISCOPLASTIC 
ANALYSIS 

1. Introduction 

As a part of the NASA Lewis’ research and development efforts for viable 
design of composite structures at elevated temperatures, a study was initiated 
to investigate efficient computational methods for the inelastic analysis of 
structural components subjected to thermal/mechanical cyclings. Our study is 
focused on large-scale computation via finite element method and heurestic nu- 
merical integration schemes for handling stiff and strongly nonlinear differ- 
ential equations associated with the constitutive modeling of high- temperature 
composites. 

It is generally known that the rate- dependent inelastic behavior (i.e. 
both plastic and creep responses) of composites can be better represented by a 
unified approach in the form of viscoplasticity theory. In this theory, the 
constitutive relations of composites are described by evolutionary equations 
which contain certain types of internal variables to model the deformation 
features such as inelastic strains, cyclic strain hardening or softening, and 
thermal ratcheting phenomenon, etc. However, the associated constitutive 
equations are known to have "stiff" regimes which present considerable numeri- 
cal difficulty in structural analysis [66]. For example, the differential 
equations, when integrated numerically, are very sensitive to the time steps 
employed. If large time steps are used in the analysis, it may result in nu- 
merical instability or inaccurate solutions [11,13,14]. If on the other hand, 
small time steps are imposed, the computational cost required for an analysis 
may become prohibitively high. In view of these considerations, a reliable 
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and efficient numerical algorithm is essential for conducting practical visco- 
plastic analysis. 

In viscoplastic analysis of structural problems via finite element meth- 
od, two levels of numerical calculations are involved: namely the structural 
and material (or local) levels. On the structural level, global equilibrium 
is enforced by using either a constant stiffness [11,13,14,67] or tangent 
stiffness approach [68-70]. On the material level, both the explicit [11,13, 
14,66,67,70] and implicit [11,69-75] integration schemes have been used to 
evaluate the viscoplastic rate equations. In addition, automatic time stepping 
in conjunction with some sort of error control may also be employed [13,14,66, 
77,78]. More recently, a uniformly- valid implicit scheme proposed by Valker 
[12,72] appears to offer great promise for obtaining stable and reliable re- 
sults in finite element analysis. 

In our first year research effort, we have conducted a numerical study on 
various integration schemes for viscoplastic materials subjected to monotonic 
as well as cyclic mechanical loadings. For this purpose, four different ma- 
terial models were implemented into a nonlinear finite element program NFAP 
[18]. These include: Valker’s model [2,19], Robinson’s model for isotropic 
materials [80], Robinson’s model for transversely isotropic materials [4], and 
Bodner’s model [81]. In addition, several numerical integration schemes have 
also been implemented into NFAP for the intended study. These include: expli- 
cit Euler forward method, explicit and implicit trapezoidal rule with or with 
out iterations, an automatic time- stepping based on the Runga-Kutta operators. 
Three numerical examples are included to demonstrate the effectiveness of the 
integration schemes considered. 
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2. Constitutive Models 


Vithin the framework of unified viscoplastic constitutive theories, vari- 
ous mathematical forms have been proposed for specific applications. In this 
section, we shall focus our attention to five particular models which appear 
to be the likely candidates for the constitutive modeling of high temperature 
composites. These are: Valker’s model, Robinson’s model for isotropic and 
transversely isotropic materials, Bodner’s model and Freed’s model. Before 
the outline of the aforementioned models, a brief description on the general 
form of the state- variable based viscoplastic equations is given below. 

General Form 

By limiting the scope of the present study to small deformation theory, 
the stress rate b in a viscoplastic material is related to its elastic strain 
rate through 


; = ?(*- j 1 ) (i) 

where b = stress rate 

e = total strain rate 
= inelastic strain rate 
D = elastic material stiffness matrix 

The inelastic strain rates are assumed to be functions of effective stresses 
in the form 



with the evolutional laws 


( 2 ) 
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( 3 ) 

* = hk If 1 ! - 7k 

where c = vector of stresses 

a = vector of back stresses 
k - drag stress 

( ) = time derivative of the quantity in question 
and h fl , hk and j a , 7k are the work- hardening and recovery functions, respec- 
tively. 

Vithin the framework of the above structure, a wide variety of viscoplas- 
tic models have been proposed. In the following, we shall list five models 
for our intended studies, 
a. Valker’s Model [79] 

The first deviatoric invariant is defined by 


I = 


\ [§ Si j - Oij] [j Si j - aij] 


1/2 


with 


a = 


2 -i. a. 

5 6 ij 


1/2 


s i j - <f\i - i ^ij (Tkk 


The inelastic strain rate is given by 


= 


T| n ^ ij " fli J ) 
l 1 


(4) 

(5) 

( 6 ) 


(7) 
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The evolution equation for the back stress rate is 




= (n 1 +n 2 )eij - (oij - o?j - n t ef j ) [(n 3 +n 4 e nsR )i+n 6 (ja k i a k i (8) 


and the drag stress rate is given by 


K = n 7 K 2 e' n7R R 


(9) 


All the material constants, i.e. n,m,ni ,n 2 ,n 3 ,n 4 ,ns ,ne ,n 7 ,a 0 and K 2 
may be temperature- dependent. K is the initial value of the drag stress and 
<r 0 indicates the amount by which the stress- strain curve is rigidly shifted 
along the stress axis. 


b. Robinson’s Model for Isotropic Materials [80] 
First, we define 


x - Sii(Si i - Qj j ) 
Vi 

(10) 


(11) 

J = 7e [(~ 23 * 8T+283S ) (777 ' t)] 

(12) 

R = Re^ 40 000 ^ ' 

(13) 

J 2 = ^ [Sij - ®ij] [Si j - aij] 

(14) 
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The inelastic strain rate is given by 



With the evolution equations being 


where: 


‘ij 


_ 2rH -I 
' UP C ij 


R gm-/?aij 


G = (G" - G°)P(x*) + G° 


f G for G > 2G° 
G " = feo for G < 2G° 


P(x) 


fO 

il^L 2 

1 


if x < - 1 
if -1 < x < 0 
if 0 < x < 1 
if x > 1 


c. Robinson’s Model for Transversely Isotropic Materials 
Assuming the existence of a potential function: 


A = ft(F,G) 


and that for small deformation the orientation direction 
plane stays the same, the two yield surfaces are defined by: 


( 15 ) 


(16) 


(17) 


(18) 

transverse 
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1 


Let 


F = & [ Xi + f Iz + 

G = KT [ X 1 + X 2 

f ( F ) = U and 


9 T 
4(4a/HI7 1 

9 

+ 4 ( 4^+17 





(19) 

(20) 


The flow law is defined by 

;?j = f (F)r t j (21) 

and the evolution law by 

a tj = h(0)eij - 7(G)Hij ( 22 ) 

where Tij are functions of the effective stress tensor Si j and the fiber 


director 

di 



F i j = F i j (^i j >didj ) 


and Hi j are given by 



Ilij - III j (fli j jdidj ) 

(23) 

with 

Si j = Si j - tti j 


d. 

Bodner’s Model [81] 


The 

inelastic strain rate has the form 



defj = /iSij 

(24) 

with 

/i 2 = D2J2 

(25) 

where 

- Z 2 (n+1) 
Dj = D 2 e n ( 3 ^ 2 ) n 
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J 2 


ySi j Si 




where 


Z - Z\ - (Zi-Zo)e 

.aVj 

m' = m 0 + mie 
Z = ffVjtdWj) - f(V°) 
dVj = Si j deij 


( 26 ) 


e. Freed’s Model [82,83] 

A viscoplastic model based on thermodynamic considerations was recently 
proposed by Freed and the model has been characterized for polycrystal- 
line metals such as aluminum, copper, nickel and tungsten. In this model 
the constitutive relations are 


S = 2/t(e - e 1 ) 


tr(<r) = 3K[tr(e) - 3a(T-T«)] 


AB = 2H j9 At 


(27) 


AD = h t At 
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where S = deviatoric stress vector 
e = deviatoric strain vector 
e* = inelastic strain vector 
(i = shear modulus 
k = bulk modulus 
B = back stress vector 
D = drag stress 
H = hardening coefficient 
The evolution equations are 


• I 


e 


= e Z 






B 

imr 


( 28 ) 


* = ll! 1 !! - \ Ilf 1 !! - * 


with a constraint equation for t or r 


7 - 1 ' \ " 2 


iz ||s|| ♦ z ||B|| a ♦ la ||B|| 
TTT 


( 29 ) 


The functional form of thermal diffusivity is adopted after Miller 


e = 


e*P(-JhO 

6xp f- Ip- + 1) 


; T > T t 
; T < T t 


( 30 ) 
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where Q = activation energy for self diffusion 
Tt = a transition temperature 
k = Boltzmann constant 

and Z is a Zener- Hollomon parameter and has the expression 

(31) 

where X = vector of effective stress 
= S - B 

At steady state, Z assumes the value 

(32) 

and A, C and n are the Garofalo creep constants. 

The limiting functions of dynamic recovery are 


J ss 


= A sinh n 


II S || 


Z = A sinh n 


II 5 

F 


L = j?D 

i = Da + (Du - Da) exp [-3(L - ||B||)] 


(33) 


Where D a , D u , and rj are material constants such that 0 < D a < D u < C. 
Further, the thermal recovery functions are given by 
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R = A 


II B ||1 


11 ? ir 

— r— J 

sinh n 

[inr] 


r = A 



sinh Q 


[II ! Ill 
~cnr 


( 34 ) 


It is noted that all the above models, except the last one, have already 
been implemented into NFAP. Freed’s model is currently being implemented 
for future studies. 


3. Finite Element Equations 

For finite element analysis, we follow the conventional incremental ap- 
proach in conjunction with an initial strain method for the treatment of vis- 
coplastic effects. In this approach, numerical calculations are segregated 
into two levels, namely, structural and local levels. On the structural 
level, the targeted solution time is divided into a sequence of (global) 
increments and the incremental structural equilibrium equations are solved for 
each time step in succession. On the local level at a material point, the 
constitutive rate equations are integrated by a subincrementing technique so 
that an incremental form is obtained. Vith this connection, the stress rate 
is related to the elastic strain rate by Eq. (1). For finite element analysis, 
it is more expedient to transform Eq. (1) into an incremental form for a typi- 
cal time step t'e[t, t + At] by writing 

= D(Ae - Ac 1 ) (35) 


where the vector of incremental stresses 
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(36) 


t+At 

= i 


a dt 


Similar expressions can be written for Ae and Ae^. Of course, efficient inte- 
gration of Eqs. (l)-(3) to obtain Eq. (35) is the main objective which will be 
discussed in the next section. 

Once the constitutive equations are defined according to Eq. (35), the 
incremental equilibrium equations for a typical displacement- based finite ele- 
ment are 


KAq = AR + AJ (37) 

in which K = the element stiffness matrix, Aq = the vector of incremental 
nodal displacements, AR = the vector of incremental applied forces, and AJ = 
the vector of incremental forces due to viscoplastic deformations, which has 
the expression 


AJ = / B T D (Ae 1 ) dv (38) 

V 

where B = the strain- displacement transformation matrix. 

In the above equation, the viscoplastic effect is treated as an equiva- 
lent (pseudo) nodal force vector appearing on the right hand side of the equi- 
librium equations. Since determination of AJ requires the knowledge of current 
state, equilibrium is approximated and is then achieved through an iterative 
process. For discussion purposes, we consider a solution step between (t n , 
and t n + At) . The equilibrium equations can be written as 


65 



At the beginning of the step: 


(39) 


KAq° = AR° 

For the i-th iteration cycle: 

KAq* = AR* + AJ 1 ; i = 1,2,... (40) 

AJi = f y B T D (Ae 1 )** 1 ^ (41) 

As shown in Eq. (9) an elastic solution, and hence the strain increment, is 
estimated at the beginning of the global time step. Based on this strain in- 
crement, the stresses, and the state variables, one can calculate the inelas- 
tic strain increment using a subincrementing procedure to be discussed later. 
Once the inelastic strain increment is obtained, the incremental pseudo force 
AJ in Eq. (41) can be evaluated. 

At this point, it is important to note that the success of equilibrium 
iterations depends upon how accurate the pseudo force AJ or the inelastic 
strain increment is calculated. Due to the "stiff" phenomenon of viscoplastic 
rate equations, the solution accuracy of Ac* is highly sensitive to the values 
of state variables evaluated at material sampling (or integration) points. 
Therefore, from the numerical standpoint, a subincrementing scheme is usually 
necessary at each integration point. If, on the other hand, the localized in- 
elastic strains are calculated incorrectly, structural equilibrium can not be 
maintained in Eq. (40). 
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4. Numerical Integration of Constitutive Rate Equations 

In this section, we present the numerical techniques for integration of 
the constitutive rate equations which are essential for the finite element 
analysis. For discussion purpose, the constitutive rate equations are 
replaced by the following expression 

y = f (y • t) (42) 

The above equation represents a system of nonlinear, first- order ordinary 
differential equations. For a multi- axial stress state, the vector y consists 
of 19 components; six Cauchy stresses <r, six back stresses a, six inelastic 
strains and one drag stress k. 

Although several numerical schemes can be used to integrate Eq. (42), 
three considerations must be taken: a) suitability for finite element imple- 
mentation, b) solution accuracy, and c) computation cost. For example, 
Gear’s multi-step methods were shown to be very effective for integration of 
viscoplastic equations under a homogeneous stress state. The methods are, 
however, not suitable for large scale finite element computing in view of two 
reasons: a) excessive computer core memory is required, and b) a special 
start-up procedure is needed. Therefore, one- step methods appear to be more 
suitable for finite element applications. 

Two classes of integration schemes are normally categorized under the 
one- step methods, that is, explicit and implicit methods. Each class has its 
advantages and drawbacks. The explicit methods are only conditionally stable, 
although less computing effort is required. On the other hand, the implicit 
methods are known to be numerically stable, but the corresponding computing 
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cost is rather high. In the following, both the explicit and implicit schemes 
and automatic stepping methods are briefly reviewed. 

4.1 Explicit Schemes 

There are several one- step explicit integration schemes available. In 
the following, we consider four different methods. 

Forward Euler Scheme : 

Ve consider a typical time interval [t n , t n + i] with the time increment 
At = t n+ i - t n . The values of y at time t n+ i may be evaluated by 

y n+ i = y n + At f n (43) 

where 

y n+ i = y(t n+ i) 


y n = y (t n ) 


(44) 


fn = f (yn, t n ) 


It is apparent that y n+ i can be calculated directly from Eq. (43) without 
involving any iterations. Eq. (43) represents the first order approximation 
to y at t n+ i • 

Euler- Cauchy’s Scheme : 

The functions of y at t n+ i are approximated by 

y = ? + j( f + f.) At (45) 

-n + 1 -'ll * -n ~ 1 
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where 


yi = y n + At • f n 


( 46 ) 


U = f(yi, t n+ i) 

Obviously, Eq. (45) represents a second-order approximation to y n + 1 . 
Trapezoidal Rule : 

By considering a second-order approximation to y n+ i» one may write 


with 


(I - £ J n ) Ay = At • f 
y n+ i = y n + Ay 


(47) 


where 




(48) 


and I = an identity matrix 

Runge-Kutta Method : 

It is a higher order integration rule. At time t n < ti < t n+ i, the "m- 
stage", k-th order Runge-Kutta method has the form 


y n + i = y n + At.^ Aj fj (yi, ti) (49) 

where Aj ’s are the weighting coefficients which have different values depend- 
ing on the order of approximation chosen; m defines the number of evaluations 
of the function f and k refers to the order of truncation error in the Taylor 
series expansion of y at t n+ i. Let E r be the truncation error. It is given 
by 
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4.2 Implicit Schemes 

Several different implicit schemes are available [71] for the solution of 
Eq. (42). Here, we will consider two different schemes: trapezoidal method 

with Newton- Raphson iterations and Valker’s integration method. 

Trapezoidal Method with Newton- Raohson Iterations : 

At time t n +i> the function y is evaluated from 


(I “ 2“ *Jn + l) Ay* +1 - yn " yn*l + (| n + + i-0,1,2. . . (51) 


At 


where Ay 1 +1 = y* + 1 - y 

- _ n+l -n+1 


and 


yn + i = y n 


(52) 


(53) 


Jn + 1 - 


'dr 

.W. 


t-tn+ l 


It is noted that the coefficient matrix of Ay i+1 on the right-hand side of Eq. 
(51) is asymmetric. If one is to use tangent stiffness method to solve the 
global equilibrium equations, cautions have to be given in the handling of the 
unsymmetric structural stiffness matrix. In general, the implicit method is 
unconditionally stable and it yields better numerical accuracy over the expli- 
cit schemes. 

Valker’s Integration Method [12,72]: 

Valker proposed a recursive integration procedure to integrate the stiff 
differential equations mentioned in the above. The relationship was obtained 



by converting the differential equations, i.e. Eq. (42), to a system of equi- 
valent integral equations. For example, we consider Eq. (7), which can be 
written as 

y + Q y = v (54) 

2 

where y = s - j ft 

s = deviatoric stress vector 
ft = back stress vector 

q = a(i-n) 

v = 2/ie - j ft 

The differential Eq. (54) can be integrated from 0 to t+At in the form 

t+At 

y(t+ 4 t) = y(t) t [<l(t+*t)-q(t)] + jr e - [q(t+it)-q(0] * d{ (55) 

t 

With an asymptotic expansion, the above equation may be approximated by 

y(t+At) = y(t) e' A ^ + ^ ^ Av(t+At) 

♦ iq < e ' iq • iq tl-e- 4 1]}[4v(t.it)-4y(t)] (56) 

where AQ = Aq(t+At) 

In Eq. (56), the evaluation of AQ and Av at t+At on the right-hand side of the 
equation requires the knowledge of y(t+At). Hence Valker’s integration method 
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is an implicit scheme. The method, according to the results presented in [12, 
72] , is very stable and highly accurate. 

4.3 Automatic Time Stepping 

In an automatic time stepping, the time increment for achieving a "spe- 
cific" solution accuracy is determined judiciously by the computer program. 
Several versions of automatic step control are available [13,66,77,78]. For 
example, a procedure based on the Runge-Kutta method was proposed in [13]. 
The key elements of an automatic step control are: a) an error estimate for 
the calculation of y at the end of each time step, and b) a criterion for 
determination of the time step At. The error involved for a particular inte- 
gration rule is given by 


U?U4t - 
St 


(57) 


where Vt+kt = f unc ti° n va l ue of y a t the time t+At evaluated by a (p+l)-th 

order integration formula 
|| || = a Euclidean norm 

A time step At in the calculation is retained if the error e computed from 
e = E/||y n || (tf-t 0 ), tf = final analysis time and t 0 = initial analysis time, 
is kept within a tolerance limit, i.e. 


e < €„ 


(58) 
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If the above criterion is not satisfied, the time step must be revised accord- 
ing to 


At' = r • At 


(59) 


where r is a non- negative coefficient, determined from 


r 


llynll 

E(tf- t 0 )_ 


1/P 


(60) 


5. Computer Implementation 

Thus far, four viscoplastic models have been implemented into NFAP for 
the intended numerical studies. These include Bodner, Valker, Robinson’s 
isotropic and transversely isotropic models. Implementation of Freed’s model 
into NFAP will be done during the second phase of this study. In addition, we 
have implemented the following numerical integration procedures into the pro- 
gram: 

- Forward Euler scheme 

- Explicit trapezoidal rule 

- Implicit trapezoidal rule 

- An automatic subincrementing control 

Valker ’s integration method will be considered in the second year study. 

In NFAP, all the viscoplastic models and the related integration proce- 
dures were programmed in a separate material module (called Model #6). Any 
future coding modifications or additions can be done locally without inter- 
fering with other parts of the program. 
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The viscoplastic material module is called by the main program during two 
distinct stages of calculations: global stiffness assemblage stage and stress 

and strain calculation stage. The material model is being driven by all the 
continuum elements: 2/D continuum, 3/D continuum, plate and shell elements. 

As shown in Fig. 21, when the material model #6 is activated, two paths of 
calculations are involved and the paths are controlled by a parameter "IND". 
That is, at the beginning of the analysis, the following subroutines are 
called for assigning workspace and array initialization: 

Subroutine Name Functionality 

INITVA Initialize VA - array which contains <r , e, a, 

and other material parameters at each integration 
point for all the elements of a particular ele- 
ment group 

M06INI Routine called by INITVA to perform the same 

function 

M06CMC Interpolate temperature dependent material con- 

stants. 

During each typical time stpe (i.e. either element stiffness calcula- 
tions, equilibrium iterations, or stress calculations), the material Model #6 
(M06STR) is called by the stress- strain driver (STSTN) and the following 
functions are performed: 

1) Formation of elastic material stiffness matrix, which may be temperature- 
dependent. 

2) Selection of a specific viscoplastic model and integration procedure. 

3) Using subincrements to calculate: A <r, Ae^, A a, etc. and update these 
quantities. 

4) Performing error checks and revising the subincrement h if necessary. 

A flow chart for the above calculations is shown in Fig. 22. 
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6. Numerical Examples 

To demonstrate the integration methods discussed in this report, three 
examples are included: 1) a uniaxial plane stress problem, 2) a simple beam 

subjected to specified displacements, and 3) a thick wall cylinder under in- 
ternal pressure. Each of these problems are described in this section. 

6.1 A Uniaxial Plane Stress Problem 

Ve consider a plane stress problem, for which a uniaxial homogeneous 
stress state is imposed using Valker (case A) and Robinson isotropic models 
(case B) , respectively. The main objective of running this problem is to test 
the convergence characteristics of various integration rules with and without 
subincrement stepping. 

Case A - Walker’s Model 

A monotonic ramp load in the form of constant strain is imposed. The 
following numerical schemes were exercised: 

a) Forward Euler Rule (Fig. 23) 

At (sec.) Number of Subincrements 

16.2 2 

16.2 4 

32.4 4 

32.4 8 

b) Automatic Subincrementing (Fig. 24) 

At = 16.2 and 32.4 sec. 
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c) Explicit Trapezoidal Rule 
At (sec.) 

16.2 

16.2 

32.2 


(Fig. 25) 

Number of Subincrements 
2 
4 
2 


d) Implicit Trapezoidal Rule (Fig. 26) 

At = 16.2, 32.4, 54 sec. 

The numerical results in the form of stress- strain plots are shown in 
Figs. 23-26. At first in Fig. 23, obviously the forward Euler’s rule is the 
least reliable unless the integration step is sufficiently small, i.e., At = 
16.2 or 32.4 sec. with number of subincrements = 4 or 8. Otherwise, numerical 
instability and/or gross error occurs in the responses. The automatic subin- 
crementing procedure with an error control (see Fig. 24) is definitely much 
more reliable regardless of the size of the global time step being used. Both 
the explicit (Fig. 25) and implicit trapezoidal rules appear to give quite 
accurate results regardless of the step size or subincrement size used. Al- 
though the computer time is somewhat more expensive as compared to the expli- 
cit schemes, the implicit method may be particularly useful for handling com- 
plex thermal- mechanical cyclings. 

Case B - Robinson’s Isotropic Model 

In this case, the material is subjected to uniaxial cyclic strains (2 1/2 
cycles as seen in Fig. 27) with three different strain rates: R = 4 x 10~ 2 , 4 
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x 10 “ 3 , 4 x 10' 4 /min. The maximum strain in all three cases was 0.327*. Both 
the forward Euler and explicit trapezoidal rules were employed. As seen in 
Figs. 28 and 29, the predictions correlate closely with the experimental data 
in [80] . It is noted however, that some numerical difficulty was experienced 
for the case of low strain rate, i.e. R = 4 x 10‘ 4 /min. This was induced by 
the discontinuous function (or shape change in the smoothing spline function) 
in Robinson’s model. Consequently, very small time steps have to be employed 
in order to reach a convergent solution. 

6.2 A Simple Beam 

A simply supported beam subjected to specified displacement at its center 
is considered. Using symmetry conditions, only one half of the beam is model- 
ed by six 8- node plane stress elements in Fig. 30. Two loading histories are 
considered in the analysis: case A - cyclic displacement with an amplitude of 
0.095 in. and loading rate of 0.1 in/min (Fig. 31), and case B - constant dis- 
placement of the same amplitude with a hold time th = 168 hrs. The maximum 
bending stress (at the center lower fiber) vs. maximum displacement of the 
beam is plotted in Fig. 32. This result is very similar to that in [80]. In 
addition, the change in maximum stress as a function of time is shown in Fig. 
33. Some stress relaxation is apparent from this plot. 

Vhen the beam is subjected to a constant displacement with a long hold 
time (th=168 hrs., significant relaxation of the maximum bending stress is 
expected as shown in Fig. 34. 
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6.3 A Thick Vail Cylinder 

A cylinder, with an internal radius 0.4 cm and external radius 0.635 cm, 
is subjected to a ramp internal pressure shown in Fig. 35. The cylinder was 
modeled by twelve 8- node axisymmetric elements under plane strain deformation. 
A global integration procedure was chosen as the following in order to achieve 
convergent solution: 

0 < t < 28 x 10* 4 hr. At = 7 x 10" 4 hr. 

28 x 10' 4 < t < 42.8 x 10-3 hr. At = 5 x 10*3 hr. 

42.8 x IQ' 3 < t < 44.3 x 10' 2 hr. At = 5 x 10' 2 hr. 


44.3 x 10-2 < t < 4.44 hr. At = 0.5 hr. 

t > 4.44 hr. At = 5 hr. 

For comparison of results, four subincrementing schemes were employed: 

Case A - automatic subincrementing with error checks 
Case B - four subincrements for all global steps 
Case C - two subincrements for all global steps 
Case D - no subincrements. 

Calculations were made from the forward Euler formula (of first order 
Runge-Kutta methods). The radial displacement at the outer wall of the cylin- 
der versus time is shown in Fig. 36 for the four cases of integration schemes. 
Apparently, the results obtained from the automatic subincrements with an er- 
ror control (limited to 0.17, according to Eq. (58) agree closely with those 
given in [80] . Solution accuracy deteriorates as the number of subincrements 
decrease. More noticeably from the result of Case D, no numerical instability 
is indicated and global equilibrium was satisfied for every incremental step. 
However, the solution is far from accurate due to the uncontrolled errors 
resulting from numerical integrations of the viscoplastic rate equations. 
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7. Conclusion 

Ve have completed the implementation of four viscoplastic models into 
NFAP for the numerical study of integration methods. Among all the integra- 
tion schemes investigated, it appears that either an automatic subincrementing 
procedure with an error control or an implicit scheme in the form of trapezoi- 
dal rule or Walker’s integral relations is more desirable. 

The major contribution of the first-year study in viscoplastic analysis 
consists of: 1) viscoplastic constitutive models and the associated integra- 
tion methods are coded as a separate module in NFAP so that any future exten- 
sion or coding modification can be done easily by a user, 2) the viscoplastic 
models can be applied to various structural components which can be represent- 
ed by 2/D continuum, 3/D continuum, or plate/shell elements, and 3) further 
understanding of the numerical schemes employed for the analysis. 

It is anticipated that our next phase of study will include: 

1) Implementation of Walker’s integration method and performing numerical 
study on various constitutive models as compared to other implicit 
schemes such as the trapezoidal rule. 

2) Implementation of Freed’s viscoplastic model for application purpose. 

3) Investigation of solution strategy for nonlinear analysis of various 
structural geometries. 

4) Development of an improved automatic stepping procedure with error con- 
trol on the basis of an implicit integration scheme. 

5) Initialization of nonlinear damage analysis. 
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PART IY: FULLY- NONLINEAR ANALYSIS CAPABILITY FOR SHELL APPLICATIONS 
1. Introduction 

The development of finite element models for the fully- nonlinear analysis 
of plates and arbitrary- curved shells, including the effects of both material 
(e.g. plasticity) and geometric nonlinearities (e.g. large displacement/ 
rotations and finite stretches), has been a very active research area over the 
years, and especially in recent years [30,31,84-89]. This extensive work was 
prompted mainly by the many conceptual, theoretical, as well as computational 
difficulties and challenging problems involved in nonlinear shell formulations 
and solutions. In addition, there is an increased industrial need to perform 
large-scale computations utilizing these general shell models, such progres- 
sive failure analysis to predict damage, buckling, post- buckling and limit 
collapse states of stiffened plates and shells, etc. 

As evidenced by numerous recent works, the degenerated- shell approach has 
apparently continued to be the most popular in such general developments. In 
particular, the important contributions here, including many variants of the 
displacement- based formulations, as well as a number of different hybrid/mixed 
developments, have been critically reviewed in [30]. However, developing a 
fully nonlinear solution procedure for arbitrarily- curved shells is a very 
demanding task. Indeed, in such a general development, it is only a minimal 
requirement to use a robust linear shell element (i.e., rank- sufficiency, 
locking free, accurate stress predictions, insensitivity to geometric dis- 
tortions, etc.) as a starting point. Additionally, several other fundamental 
issues must be carefully investigated. Confining attention to the quasistatic 
problem, the following are of utmost importance from a theoretical standpoint: 
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(i) consistent linearization of the underlying weak, or variational, form of 
the governing equations [88-90]; (ii) treatment of large rotations , in both 
stiffness derivation as well as the configuration update procedure [31,91,92]; 
(iii) use of objective measures of stress and strain, and their rates, suit- 
able for the particular form of the nonlinear material model employed (e.g. 
plasticity, viscoplasticity, etc.) [93-98]; and (iv) use of a valid method 
for stress integration (update) that maintains incremental objectivity in the 
presence of large rotations and "non- inf initesimal" stretches, when dealing 
with rate- type constitutive equations [75-99]. 

The objective in the present paper is to briefly review some of our 
recent work dealing with the formulation of a general theoretical framework 
for the analysis of finitely stretched and rotated shells using our hybrid/ 
mixed method [30] . In particular, this is utilized here as a basis for 
deriving one representative nonlinear model of this type; i.e., HMSH5 (four- 
noded quadrilateral) . Numerical results illustrating the performance of the 
element in some test cases are also reported in order to demonstrate the 
effectiveness of various proposed solution schemes. More examples can be 
found in [30] . 

Ve finally note that the present mixed element HMSH5 is "extended" 
nonlinear version of its linear counterparts; e.g. see [27] for applications 
to isotropic shells, and Part II for extensions to laminated plates/shells. 
The crucial point in its mixed formulation concerns the judicious selection 
of the parameters in the polynomial strain approximation. In particular, this 
was facilitated by the use of a set of bubble functions (i.e., kinematic 
degrees of freedom associated with an interior fifth node) in HMSI5. 
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For conciseness, the scope of the present discussion is limited to quasi- 
static problems. In addition, special instability phenomena and post- buckling 
responses are not addressed. Further, in all the applications presented 
later, we utilize a "semi- linear" elastic isotropic model, thus enabling the 
comparison of our results with other independent solutions available in the 
literature. 

2. Geometric and Kinematic Descriptions 
2.1 Basic Hypotheses 

The main assumptions underlying the present mixed- element developments 
are summarized as follows: 

(i) straight and inextensible shell "normals" (or "thickness" fibers). 

(ii) zero "through- thickness" normal stress (i.e., plane- stress assumption). 

(iii) "small" transverse shear strains, but large membrane and bending strain 
components. 

The first two simplifying hypotheses (i) and (ii) are typical in any de- 
generated shell model [e.g., 84], based on classical Mindlin/Reissner plate 
theories. However, note that, even though assumption (i) is utilized as a 
basis for deriving the element governing "stiffness" equations , fiber "thin- 
ning" effects due to large membrane strains can still be recovered in the 
computations; e.g., using the incremental "average- thickness" update procedure 
in [86] . 

On the other hand, certain important implications resulting from the ad- 
ditional assumption (iii) can be exploited to develop simplified and effective 
numerical algorithms [30]. From a practical standpoint, its justification is 
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that "large” transverse shear strains are not typical in common thin or moder- 
ately- thick composite shell structures. 

2.2 Coordinate Reference Frames 

Three different types of Cartesian reference frames are defined here for 
convenience in subsequent derivations: 

(1) The fixed or global reference frame, with its orthonormal base vectors 
e^(i=l,2,3), is used to define element geometry and its translational 
displacement DOF (degrees of freedom) . 

(2) A unique local fiber coordinate system is constructed at each node, with 
associated base vectors e|(i=l,2,3), where e^ coincides with the nodal 
fiber. During the incremental analysis, this orthogonal fiber triad is 
continuously updated and used as a moving basis for defining the "finite" 
rotational DOF at the node. 

(3) A local lamina system at each integration point in an element, 

i l 

e^(i=l,2,3), with e^ taken to be normal to the lamina surface, and the 

/ / 

other two in- plane lamina- tangent vectors e^ and selected as in [30] 
(see Fig. 37) . 

This latter lamina basis rotates rigidly as the element deforms, and it 

is found to be most convenient for invoking the plane stress assumption in the 
i 

2 * direction in all configurations as well as for ensuring the invariant pro- 
perty in the interpolation of the assumed strain field given later. 


83 


2.3. Geometry and Kinematics 

Vith its midsurface taken as the reference surface, the position vector 
to an arbitrary point in the shell element at any time instant "t" is defined 
in terms of the natural coordinates (r,s,t) as follows: 


t x = S 6 N k t x k + l S 6 t h,_ N r VW 


k=l 


k=l 


‘k k -3 


( 2 . 1 ) 


where ^x k , and ^h k are the position vector, components of a unit 

pseudonormal (fiber) vector, and fiber- dimension (shell- thickness) parameter, 
respectively, at nodal point k on the reference surface. The N k (r,s) are the 
two-dimensional shape functions associated with node k [27], and n the number 
of nodes per element. 

There are five DOF per node used to parameterize the element configura- 
tion in the shell space: three global translations (u,v,w), and two "fiber" 

f f 

rotations and e ^ about axes e^ and Note that, following [30], these 
latter nodal rotation freedoms are appropriately viewed here as generalized 
finite rotational coordinates of the Rodrigues- Euler type [91,92], thus pro- 
viding a convenient, singularity- free (no drilling freedoms) parameterization 
for the director fields of the elements. In total, an HMSH5 element thus con- 
tains 25 DOF. 

In the context of the present incremental analysis, three successive con- 
figurations, at time "0" (initial), "t" (current) and "(t+At)" (incremented or 
neighboring) are considered. Then, the total and incremental displacement 
fields of the element can be expressed as 
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( 2 . 2 ) 



(2.3) 


The relation in Eq. (2.2) is directly employed to evaluate the total ele- 
ment displacements and their derivatives (i.e., total geometric "Almansi" 
strains in Eq. 3.1). In this, components of the "updated" director vectors, 
e^, in the configuration "t" are determined in terms of nodal rotations using 
the (geometrically- exact) update procedure of Sec. 4.4. 

On the other hand, the relation in Eq. (2.3) provides the basis for de- 
riving the "linearized" governing equations of the element. To this end, the 
second term must first be expanded in terms of nodal rotation increments, and 
here we utilize the following "linearized" kinematic approximation: 

n e l n e 

Au = E N k Au k + 5 E \ N k (-Ae k t e^ k ^ + Ae t e^ k ^) (2.4) 

for the stiffness evaluation of both elements considered. Ve particularly 
note that this leads to the well-known expression (i.e., similar to the "true" 
continuum case [84]) of a single geometric stiffness matrix (see Sec. 4.2). 


n r n 

A v G w a v e t. « /t+Atf(k) t f(k)\ 

' S 1 N k A -k + 2 J S 1 h k N k ( ?3 ?3 ) 

k=l k=l 


?k * S ■ \ and 4 ?k = t,At ?k - ‘?k 
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Remark 2,1 Alternative forms for nonlinear kinematic approximations of the 
rotational terra in Eq. (2.3) have been also employed in several recent studies 
[e.g. 88,89]. This led to the introduction of various additional geometric 
stiffness contributions, which were found to be necessary in order to attain 
quadratic rate of (asymptotic) convergence in Newton- Raphson iterative schemes 
for solutions of large- rotation shell problems. However, as was demonstrated 
in [30], even with the single geometric stiffness based on (2.4) above, the 
same "good" convergence rate is also exhibited by the present mixed model 
HMSH5; the crucial point here is the "exact" updating of the orientations of 
nodal fiber triads (Sec. 4.4). 

Remark 2.2 In addition to its computational efficiency, it was also proved in 
[30] that the geometrix stiffness of HMSH5 possesses a "second- order" accur- 
acy, thus making the element capable of appropriate modeling of instability 
modes; e.g. creep buckling. 


3. Variational Principle 

A modified Hellinger- Reissner variational principle [e.g. 30,100] pro- 
vides the starting point for the present incremental, step-by-step analysis. 
This takes the following form in updated Lagrangian (UL) description with "t" 
as reference; i.e., ^Ar^ = 0, 


Ar, 


HR 


Afp W A TO A /p A 

- 2 Ae A cAe + <r Ae + Ae cAe - Ae A c(e-e) 


dV - AV 


(3.1) 
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where Ae and Ae = Ae (linear) + Aj/ (nonlinear) are, respectively, independent- 
ly-assumed and geometric (from displacements) Green strain increments; e and e 
the corresponding total Alamansi strains; c the material stiffness; A<r = cAe 
the Truesdell stress increment; o the true (Cauchy) stress; and AV is the work 
of prescribed forces. The last term in the bracket of Eq. (3.1) is due to 
compatibility- mismatch [5,45]. 

Note that, for the purpose of implementing element HMSH5, all strain/ 
stress vectors in the above will be defined with respect to the lamina basis 
at "t". In particular, this implies the use of a "reduced" (5x5) material 
matrix, in accordance with assumption (ii) of Sec. 2.1. 

4. Finite Element Formulation 
4.1 Strain- Field Discretization 

In addition to displacement interpolation, a polynomial for Ae in the 
present mixed element is also needed. To this end, we utilize the same spe- 
cific "least- order" polynomial strain approximations for these elements pro- 
posed previously for linear analysis [30]. Thus, in tensor- component form, 
the incremental lamina strains Ae' for undistorted geometry are, for HMSH5 
element: 


Ae n 

tl 

+ 

/?2 r + 

/? 3 s + f(/? 4 + 

4> r * 

Ae 22 

■V 

P 8 r + 

Pg s + HP 10 

+ ^ll r + ^12 S ) 

Ae 12 

= ^13 

+ P 14 l 



Ae 23 

^15 

+ ^16 r 

* 


Ae 13 

= ^18 

+ ^19 S 

* V 
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Note that the above strain distribution needs to be further modified to 
account for geometric distortions; namely the "important" in- plane lamina 
(skewness) distortion. To this end, we use a "constant" Jacobian transforma- 
tion as described in [25,30]. Vith this, we can finally write 

Ae = P A£ (4.2) 

where A/? are generalized strain parameters, and P are "modified" strain- 
interpolation functions. 

4.2 Element Stiffness Equations 

In view of Eqs. (2.4) and (4.2), the appropriate "linearized" form of the 
variational principle in (3.1) is simply obtained as in the convectional 
"true" continuum case. This yields, after invoking the stationarity condi- 

tions with respect to A/?, and then Aq, 

A/? = f 1 [G Aq + (Z x - Z 2 )] (4.3) 

where 

H = / P T c P dV ; G = / P T c B T dV (4.4) 

V ' ' " • V • • - L 

Z t = J P T c e dV ; Z 2 = | P T c e dV (4.5) 
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as well as the desired final stiffness relationships: 

(?L + *nl) a ? * 9 • (9i ♦ y 


(4.6a) 


?L - f ! 1 9 ! ?NL = J ?NL i ?lfl dV < 4 -« b ) 

= J ?L ? dV ; q 2 = f 5 1 (?1 ■ ?2> (4.6c) 


See details in [27]. The is the element "linear" (or material) stiffness, 
its geometric stiffness (shown in [30] to exhibit second- order "accur- 
acy"), and the right-hand side of Eq. (4.6a) includes "correction" terms due 
to both equilibrium imbalance (Q- Q^) as well as compatibility mismatch terms 
in q 2 . 


4.3 Solution Procedure 

Once assembled, the linearized equations above are utilized in the fol- 
lowing global incremental- iterative full Newton- Raphson scheme: 

t t+4t K Nt )( n ) 4q(" n ) = UAt Q - t+At (Q 1 ♦ Q 2 )W (4.7) 

for the solution for the incremental nodal displacements in the (n+l)th iter- 
ation within the time step "t"-*"t+At". A displacement- type convergence cri- 
terion is adopted here (0.001 tolerance for the ratio of norms of incremental 
to total nodal displacement vectors). 
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Remark 4.1 In applications to viscoplastic analyses, the above variable 
stiffness solution scheme calls for the formulation of a (material) constitu- 
tive tangent stiffness matrix C (Eqs. 4. 6,4. 7). Although such matrices are 
not readily available at present for "general" unified viscoplastic models 
[e.g. 4,82,83], several attempts have been made [68-70] in the context of 
other viscoplasticity theories. This aspect will be further investigated in 
our future work. 

Remark 4.2 Instead, the more conventional initial strain or constant stiff- 
ness , viscoplastic solution approach is obtained here by simply utilizing an 
elastic C in all Eqs. (4.3)- (4.6) above, together with an additional "load" 
term (i.e., Q 3 - vector) in the parenthetical term on the right-hand side of Eq. 
(4.6a) to account for inelastic strain- rate effects. 


4.4 Large- Rotation Configuration Update 

Configuration update involves the calculations of (i) new nodal coordi- 
nates, as well as (ii) orientations of the associated fiber triads. Although 
(i) is trivial, (ii) is complicated by the non- vectorial character of finite 
space rotations. Here, we make use of the so-called exponential mapping 
algorithm for rotational transformation of vectors [30,90-92]; i.e. at the end 
of the (n+ 1 )"^ iteration, 


e f ( k ) 
ri 

+ 

» 

e f 00 

-2 

II 

e f ( k ) 
?3 J 
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(4.9) 


1 cos ||Ae|| sin ||Aejj 

Ml' 2 ; 62 = Mil 


||A ? || « (« 2 ♦ 


where the local components of the rotation pseudovector Ae^ n+ ^ along the 
fiber axes e^) ( a t configuration "n") are conveniently defined as (n,/?,0). 


4.5 Strain Update 

The calculation of the updated independent (Almansi) strain field e at 
the quadrature points in the element is needed to evaluate the iterative com- 
patibility mismatch "force" vector Qg. A simple "approximate" procedure [30] 
is utilized here. In this, the updated e is obtained by a push- forward 
transformation [90] for the (covariant) tensor components of the total (in- 
cremented) strains using the relative deformation gradient for configurations 
(n+1) and (n) . For convenience, the latter is written as follows, using the 
polar- decomposition theorem [90], 








(4.10) 


where R and U are rigid- rotation and pure stretch tensors, respectively. With 
appropriate reference to lamina coordinate systems in different configura- 
tions, we may then show that 

? (n * 1) = f„ii(? (n) * A ? (n+1) ) ?»!i ( 4 - u ) 
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where now all components of tensors on the right-hand side are referred to 
lamina coordinates in (n), whereas e^ n+1 ^ are taken with reference to lamina 
axes in the new configuration (n+1). 

Remark 4.3 In keeping with the present mixed formulation, it is crucial to 
determine the relative stretches (as well as all other pure "strain- like" 
quantities) from the "true" strains. To this end, we use the following ap- 
proximation [30]: 


0 n+1 = (1 * |l|)I*(l - lx|)4e (n * 1) - j(l-lf)Ae(” +1 )Ae(" +1 ) (4.12) 

where l| = (i=l,2,3) are the invariants of Ae( n+ ^. 

This formula obviates the need for "costly" procedures to formally obtain 
the square- root of a positive- definite matrix, while it still maintains the 
condition of same principle axes for and Ae( n+1 ). 

Remark 4 . 4 However, when needed later (Sec. 5.1), must be calculated 
from element displacement field. For "relatively" small strain increments, it 
can be approximated [30] by the coordinate- transformation matrix for lamina 
bases in configurations (n) and (n+1). 
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5. Stress Update 

For a class of moderately large- strain constitutive models, a general spatial 
rate- form, consistent with the variational statement in Eq. (3.1), is consi- 
dered here (a superposed dot indicates a material time derivative) 


O 

a 


ij 


c ijkl d kl 



{ i - f + (d kk ) 


(5.1) 


where a is the (objective) Truesdell rate of Cauchy stress a, d the (spatial) 
deformation- rate tensor, l the velocity gradient, and c may generally depend 
on stress and/or deformation history (e.g. plasticity). 


Remark 5.1 Considering the range of "non- inf initesimal" strains anticipated 
in viscoplastic/damage studies for composite space- engine components (i.e., 
2- 5%) , the precise distinction between various "objective" stress rates be- 
comes unimportant from the standpoint of elastic modeling of material behavi- 
or. However, from the viewpoints of both numerical implementation as well as 
solution accuracy/convergence, the use of Truesdell rate in Eq. (5.1) is of 
great advantage in the present mixed formulation, i.e. leading to "efficient" 
integration algorithm in Eq. (5.4) and "improved" numerical convergence (e.g. 
Sec. 6.3). 

5.1 The Basic Stress- Integration Algorithm 

In addition to the requirements of numerical stability and consistency , a 
time- stepping scheme for stress updating with finite strains/rotations must 
also satisfy the important condition of the incremental objectivity, [31,95]. 


C -3~ 
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In particular, because of its desirable numerical characteristics, the gener- 
alized mid- point (trapezoidal) rule has been utilized extensively to formulate 
"objective" stress- updating procedures. In this section, the method is adapt- 
ed for the present mixed shell elements. To this end, the following defini- 
tions are introduced: 


F = R U , J = det F 
-n+a -n+a -n+a ’ n+a -n+a 


(5.2) 


4 ?<”* 1 / 2 ) ‘ 4t ?n+l/2 ‘ -n+1/2 ?n+l/2 > ?n+l/2 “ 4t ^n+1/2 ( 5-3 ) 


in which c n+ jy2 is the fourth- order material moduli tensor at the mid- point 
configuration t n+1 ^ 2 > At = (t n+ j- t n ) , and 0 < a < 1 (with F fl = I, = 1, 
where I is the unit tensor) . 

Vith the Eqs. (5.2) and (5.3) the stress update algorithm takes the 
following general form, giving now, the updated stress components a^ n+ ^, 
directly referred to the updated t lamina basis [30]: 


-(n+1) . j-1 r («r( n ) + J 
~ J n+l-n+l'- J 


,-l 


4 ; (n+1/2) r„!i/2)!„ + i ( 5 - 4 > 


'n+1/2 ?n+l/2 ' -n+l/2'-n+l 


r ( n+1 ) 


in which all tensor components, except <n ‘ , are measured in the t - lamina 
basis. 


Remark 5.2 For the important common case of infinitesimal 

rotation analysis, F„ ~ I and J ~ 1 (0 < a < 1), and Eq. 

J 7 -n+a - n+a v ^ 

reduces to (<r n+ * in lamina basis at t ,): 
v - n+l 7 


strain/ large 
(5.4) simply 
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(5.5) 


,(" +1 ) = »(») ♦ A> +1 ) I A> +1 ) = c Ae( n+1 ) 


5.2 Approximations 

Now, with Eq. (4.12) available, our remaining task is to determine dis- 
crete approximations for F n+1 ^ 2 and D n+1 / 2 from known (n)- and (n+1)- configur- 
ations, and Ae( n+1 ). To this end, certain additional assumptions are needed 
concerning the deformation ("integration") path during the incremental step. 
Here, two common assumptions are utilized: (i) straightline deformation path 
and constant velocity for each material point, or (ii) constant material rota- 
tion and fiber extension rates. 


For Assumption (i), we can easily show that 


?n+l/2 = 5<?n*l + !> i ?n+l/2 = 2 [<?n*l ‘ P(?„*l * P' 1 ]* < 5 ' 6 ) 

in which "s" indicates the symmetric part, and again we emphasize that F fl+ ^ is 
calculated from its component tensors R n+ ^ (Remark 4.4) and D n+ ^ (Eq. 4.12). 

For Assumption (ii) here, after "somewhat" lengthy derivation, we find: 


where 


■ C tilt t ?„ +1/2 - Krt/jtfc ! n+ i]«Li/2 


An = If [An A]« T ; £ n+1 = M A M T 


(5.7) 

(5.8) 
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The M is an orthogonal matrix containing eigenvectors of Sn+1’ and diagonal 
matrices A and [£n A] contain its eigenvalues (principal stretches) (i= 
1,2,3), and their natural logarithms, respectively. 

Remarks 5.3 Instead of Eq. (5.8), the following first-order, Pade-type, ap- 
proximation may be utilized: 

!„*1 = 2 (! n + l - !)(!».!* I )' 1 ( 5 - 9 ) 

6. Sample Applications 

Ve consider here three numerical simulations for HMSH5. All results are 
obtained assuming isotropic semi- linear elastic material behavior, and except 
for some large- strain applications in Sec. 6.3, all other problems were solved 
using the small- strain/large- rotation assumption. Ve refer to [30] for addi- 
tional examples. 

6.1 Clamped Square Plate Under Uniform Load 

This problem is taken from [31] . Figure 38 depicts the results of HMSE5 
utilizing the same 4x4 mesh as for the 4-noded (reduced- integration) model 
Q4- UI in [31]. It is interesting to note that the average number of itera- 
tions per load step (about 3) for the two elements was the same. 

6.2 A Pinched Cylinder 

This problem has often been used in establishing the viability of shell 
elements in linear analysis [27]. But to our knowledge, there is currently no 
"benchmark" solution available for its nonlinear response. Only one- eighth of 
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the shell was modeled (because of symmetry) using a 16x16 element mesh for 
HMSH5 (Fig. 39). Ten equal load increments were applied to arrive at the full 
load of 750 (about four times the load utilized in linear solution). A total 
of 39 iterations were needed for HMSH5. 

6.3 A Pinched Hemisphere 

This is another "obstacle test" in linear analysis. It is also an excel- 
lent test of the ability of an element to handle "truly" three-dimensional 
finite rotations. Using symmetry, one quadrant of the shell is modeled with 
10x10 mesh for HMSH5. The total load F=100 is first considered here for the 
small/large rotation analysis (recall that F=1 is typically used for the lin- 
ear case) . Only two loading steps are required for HMSH5 from F=0-*10 and 
F=10-*100 with a total of 13 iterations. Results are shown in Fig. 40, to- 
gether with the solution reported in [89] using the resultant- based shell 
element, Q4S , with a 16x16 mesh. 

In addition, a large strain analysis using the algorithm of Eqs. (5.4) 
and (5.6) and the same 10x10 mesh for HMSH5 was also performed for the above 
shell subjected to a very high load level F=90O. Ve utilized two and 8 equal 
load steps for F=0-*100 and F=100->900, respectively. Significant deformations 
occurred in this case; e.g. at the final load level (F=900) , u^=4.922 and 
Vg=12.709, with a total rotation of nearly 123 degrees at B. A summary of the 
results is given in Table 11, where we also report the results obtained from 
an independent inf initesimal- strain solution version for comparison. Note 
that the total (accumulated) strains at F=900 are now of the order of 37.. 
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As evident from Table 11, although there are no discernible differences 
between the two solutions above for displacements and rotations at most load 
levels, the "large- strain" update procedure produced significant reductions in 
the number of equilibrium iterations required for convergence at higher load 
(strain) levels; e.g. from 23 to 16 iterations for F=500-+900. This is, of 
course, due to the more accurate stress calculated in this latter case, and 
will become even more important in cases when "stress- dependent" viscoplastic 
models are utilized. 

7. Conclusions 

Ve presented here the fully- nonlinear formulation of curved shells using the 
simple mixed model HMSH5. Several noteworthy aspects are included. A careful 
selection is made for the polynomial functions in the strain assumption, thus 
leading to robust elements (i.e., kinematically stable, free from locking, 
etc.). A geometrically- exact procedure is utilized for element configuration 
update with finite nodal rotations. Even with a "single" geometric stiffness 
matrix, this was shown to be capable of attaining quadratic convergence rate 
in practical shell applications [30] . It is also the complete geometric ma- 
trix needed for adequate modeling of "significant" instability modes; e.g. in 
creep buckling studies of composites under thermomechanical loadings. 

For updating of the spatial stress field in the presence of non- infini- 
tesimal strains, "objective" generalized midpoint schemes were developed by 
making use of the polar- decomposition of the deformation gradient. These are 
in keeping with the underlying mixed method. 

Several numerical simulations have been presented to demonstrate the ef- 
fectiveness and practical usefulness of the proposed formulation. 
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Our future work will focus on extensive testing of the validity of the 
proposed solution algorithms in the context of life prediction studies of 
high- temperature composite plates and shells. In particular, a variable 
stiffness formulation for unified viscoplastic models with damage will be 
developed and used in several "benchmark" solutions of static and dynamic 
problems. 
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Table 1 


A Cantilever Laminated Beam with In- Plane Load 


Laminate 

Deflection 

at 

TiD 

Exact 

■UJ.UtkiBAJlM 

TRIPLT 


r 90] s 

II 


0.0864 

1.0956 

0.3543 

0.2353 

0.1673 


Table 2 

Mesh Convergence Study for Clamped 
1 Layer Rectangular Plate 
Normalized Center Deflection 


Laminate 


Present 

F.E. 


TRIPLT 
(8 X-8)_ 

(4 x 4) 

(6 x 6) 

KEm 


[5] 

10.4598 

10.9851 

10.9592 

10.8879 

10.5375 

L J 

(-0.747.) 

(+4.257,) 

(+4.07.) 

(+3.337.) 


[15] 

9.3339 

9.9972 

10.0391 

10.0215 

9.4455 

(-1.187.) 

(+5.847.) 

(+6.287.) 

(+6.107.) 


[25] 

5.9319 

6.5407 

6.5245 

6.5132 

6.0856 

L J 

(-2.537.) 

(+7.487.) 

(+7.217.) 

(+7.037.) 


[35] 

2.9553 

2.9991 

2.9269 

2.9226 

2.8985 


(+1.967.) 

(+3.477.) 

(+0.987.) 

(+0.837.) 


[45] 

1.6414 

1.3695 

1.3644 

1.3820 

1.4139 


(+16.097.) 

(-3.147.) 

(-3.507.) 

(-2.267.) 


[75] 

1.1334 

0.7911 

0.9000 

0.9133 

0.9134 


(+24.097.) 

(- 13.397.) 

(-1.467.) 

(0.7.) 


[90] 

0.9962 

0.6608 

0.8116 

0.8073 

0.8017 


■Z HEsEM 

(- 17.577.) 

(+1.237.) _ 

(+0.707.) 

A K 

a a r 


D.O.F. 125 245 405 605 1215 
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Table 3 


Clamped Rectangular Plate 
with Uniform Pressure 



Normalized 


Center 

Laminate 

Def] 

Lection, w 


TRIPLT 

Present F.E. 

■HI 


10.537 

10.888 


15 


9.4455 

10.022 


m 


6.0856 

6.5132 


45 


2.8985 

2.9226 


60 


1.4139 

1.3820 


75 


0.9134 

0.9133 


Elil 


0.8017 

0.8073 


Table 4 

Clamped Square Plate 
with Uniform Pressure 


Laminate 

Normalized 

Center 

Def lection, w 

Exact 

Present F.E. 

M0H3T 


05 


0.0946 

0.1040 

0.1083 


[±25‘ 


0.2355 

0.2602 

0.2572 


±35 


0.2763 

0.2914 

0.2844 


±45 


0.2890 

0.3013 

0.2929 


Table 5 

Simply Supported Square Plate 
with Uniform Pressure 


Laminate 

Normalized 

Center 

Deflection, w 

Exact 

Present F.E. 

TRIPLT 

V2R 1 


05 





mam 


[±25" 





1 


±35 


0.945 

1.041 

0.952 



±45 


0.9150 

1.033 

0.922 
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Table 6 


Bending Moment M for a Simply Supported 

AA 

2 Layer Square Plate 


Laminate 

M 

,xx 

Exact 


TRIPLT 

V2R 

Wm 

jKj 


1318. 

1509. 

1327. 

1396.2 

■ 

fUl 


1142. 

1137.6 


1176.8 

1 


1 

843.6 

842.5 

851.3 

855.6 

I 


II 

564.6 

591.9 

573.1 

567.8 

■ 

7ijt 

|| 

368.1 

422.17 

375.3 

371.3 


Table 7 

Bending Moment M yy for a Simply Supported 
2 Layer Square Plate 


Laminate 

M 

vy 

Exact 

■I&SS399H 

TRIPLT 

V2R 


[±5] 


34.2 

37.56 

34.47 

34.99 


[±15‘ 


123.4 

123.64 

124.4 

126.4 


'±25' 



221.37 

228.4 

221.5 


±35 



314.02 




±45 


368.1 

422.17 

375.3 

371.3 1 


no 


















Table 8 


Mesh Convergence Study for Isotropic 
Cantilevered Plate 


Freauencv 


Present 

F.E. 



(4 x 4) 

(6x6) | 

(8 x 8) 

liEBSEIB 

Ref. f65l 

1 

3.34 

3.39 

3.41 

3.42 

|^KKTM| 

2 

13.32 

13.82 




3 

19.46 

20.43 


20.98 


4 

40.94 

43.58 

44.91 


: JiisRSi5S 

5 

51.22 

55.32 

56.98 

57.82 

60.50 

6 

71.62 

78.81 

83.08 

85.38 

92.30 


Table 9 

Anisotropic Cantilever Beam Solution 
Fundamental Frequency, u\ 


Fiber 

Oden el 

’H'wni— 

Present 

Direction 

ill 

(id: 

F.E. 


14.351 

14.130 

14.232 



13.020 

13.086 

45 

12.873 

12.679 

12.757 

60 


13.020 

13.086 

90 

14.351 

14.130 

14.232 


Table 10 

Anisotropic Cantilever Beam Solution 
Second Frequency, ui 


Fiber 


i.al [24] 

Present 


(i) 

an 

F.E. 

0 

85.950 

83.675 

84.484 

30 

80.673 

78.264 

78.147 

45 

78.379 

76.499 

76.389 

60 

80.379 

78.264 

78.147 

90 

85.950 

83.675 

84.484 


ill 





















Table 11 


Comparison of Results for the Solutions of the 
Pinched Hemisphere (10 x 10 mesh) ; see Fig. 40 


Load 

Infinitesimal- Strain Solution 1 

Laree- Strain Solution 

F 

U A 

V B 

n e 4 ii 

ll®B« 

Iterations 

U A 

V B 



Iterations 




(rad) 

(rad) 






50 

2.562 

3.797 



6 

2.562 

3.797 

0.543 


6 

100 

3.295 

5.713 



5 

3.295 

5.711 

0.727 


5 

200 

3.867 

7.637 


1.470 

5 

3.867 

7.636 

0.907 

1.469 

5 

300 

4.162 

8.721 

1.023 

1.663 

4 

4.164 

8.725 

1.023 


4 

400 

4.426 

9.712 

1.138 

1.802 

5 

4.429 

9.721 

1.139 

1.803 

5 

500 

4.650 

10.763 

1.234 

1.920 

6 

4.652 

10.772 

1.235 

1.921 

6 

600 

4.782 

11.559 

1.269 

2.008 

6 

4.782 

11.563 

1.270 

2.009 

4 

700 

4.851 

12.070 

1.260 

2.069 

5 

4.851 

12.067 

1.260 

2.069 

4 

800 

4.894 

12.432 

1.224 

2.114 

5 

4.893 

12.422 

1.226 

2.114 

4 

900 

4.925 

12.739 

1.146 

2.153 

7 

4.922 

12.709 

1.165 

2.150 
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Figure 1: 4- Node element at mid- plane of multi-layered composite 
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Figure 2: Local fiber orientation with respect to element coordinates 
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Figure 3: Multi-layer composite 







Figure 5: 
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Clamped rectangular plate 
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Figure 9: [± 15] <r x stress distribution 
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Figure 12: r xz after strain re-distribution, approach 1 
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Figure 13: r xz distribution case 1, approach 2 [0/90/0] laminate 
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STRESS 


Figure 14: r xz distribution case 2, approach 2 [0/90/0] laminate 


124 





3.00 H.OO 5.00 6.00 


STRESS 


e 3, approach 2 [0/90/0] laminate 


125 




LOCfiT ION 



SHEAR STRESS 


Figure 16: r zz distribution approach 3 [0/90/0] laminate 
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Figure 18: Mode shapes for first and second frequencies 
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Figure 19: Effect of thickness on the fundamental frequency for 

two and eight- layer [0/90/...] laminates 
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20. 



ndamental frequencies vs. plate aspect ratio for 
/90/0] laminates 
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Figure 20 : Fui 
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Figure 202: Fundamental frequencies vs. plate aspect ratio for 
four- layer laminates 
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Main Driver of Viscop: 
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Fig. 23. Uniaxial Stress-Strain Response of Walker's Model 

Using Different Step Sizes 
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Uniaxial Stress-Strain Response of Walker's Model 
Using Self-Adaptive (Automatic) Stepping 
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Uniaxial Stress-Strain Response of Walker's Model 
Using Explicit Trapezoidal Rule 
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Fig. 26. Uniaxial Stress-Strain Response of Walker's Model 

Using Implicit Trapezoidal Rule 
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tress-Strain Resi 
Using a Forwar 











Fig. 29. Comparison of Cyclic Stress-Strain Responses of 

Robinson's Model Using Euler and Trapezoidal Rules 
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Fig. 30. Finite Element Model of a Simple Beam 
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Cyclic Displacement Loading for a Simple Beam 





PLflNC 3TK 



143 





rLflNC 3TRC35 BNHLr3l3 
8 CONST. SUBIN. 



144 


Maximum Stress History in the Simply Supported Beam 
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Stress Relaxation in a Beam Under Constant Displacement 
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36. Radial Displacement at the Outer Wall of 
the Cylinder with Different Incrementing 
Schemes 
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Fig. 38. A Clamped Square Plate Under Uniform Load 
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Fig. 40. A Pinched Hemispherical Shell 
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